Input

Analysis for FDCR manuscript describing the temporal dynamics of drug cue reactivity with replication. ## Read Data

#function to be able to read in an .RData file and assign it to a variable
loadRData <- function(fileName){
  #loads an RData file, and returns it
  load(fileName)
  get(ls()[ls() != "fileName"])
}
idps_neurocaps <- loadRData('../paper-dcr-temporaldynamics/data/neurocaps/idps-neurocaps-2-13-2020.RData')
idps_neurocaps$motion <- idps_neurocaps$dcr_motaveraw_r1
#remove subjects due to excessive motion
idps_neurocaps <- idps_neurocaps[idps_neurocaps$motion < 0.3 & (!idps_neurocaps$neurocaps_exclude),]
idps_neurocaps_ocr <- idps_neurocaps[idps_neurocaps$dcr_version == 'ocr',]
idps_neurocaps_mcr <- idps_neurocaps[idps_neurocaps$dcr_version == 'mcr',]
idps_tacs <- loadRData('../paper-dcr-temporaldynamics/data/tacs/idps-tacs-2-13-2020.RData')
idps_tacs$motion <- idps_tacs$dcr_tpre_motaveraw_r1
idps_tacs <- idps_tacs[idps_tacs$motion < 0.3 & (!idps_tacs$tacs_exclude),]
idps_tdcs <- loadRData('../paper-dcr-temporaldynamics/data/tdcs/idps-tdcs-2-13-2020.RData')
idps_tdcs$motion <- idps_tdcs$dcr_tpre_motaveraw_r1
idps_tdcs <- idps_tdcs[idps_tdcs$motion < 0.3 & (!idps_tdcs$tdcs_exclude),]

#only use neurocaps OCR subjects who also did tACS
idps_neurocaps_ocr <- idps_neurocaps_ocr[idps_neurocaps_ocr$id %in% idps_tacs$id,]

#only use tacs subjects with good neurocaps data
idps_tacs <- idps_tacs[idps_tacs$id %in% idps_neurocaps_ocr$id,]

#idps_tacs$dcr_tpre_response_craving_[0-7]
#idp_tdcs
#dcr_tpre_response_box_rt_0
#dcr_tpre_response_craving_rt_0
#idps_neurocaps_ocr$dcr_response_box_craving[0-7]


n_tdcs <- length(unique(idps_tdcs$id))
n_neurocaps_mcr <- length(unique(idps_neurocaps_mcr$id))
n_neurocaps_ocr <- length(unique(idps_neurocaps_ocr$id))
n_tacs <- length(unique(idps_tacs$id))


print('Subjects In Each Dataset')
## [1] "Subjects In Each Dataset"
print('tDCS MCR')
## [1] "tDCS MCR"
print(n_tdcs)
## [1] 65
print('Neurocaps MCR')
## [1] "Neurocaps MCR"
print(n_neurocaps_mcr)
## [1] 29
print('Neurocaps OCR')
## [1] "Neurocaps OCR"
print(n_neurocaps_ocr)
## [1] 22
print('tACS OCR')
## [1] "tACS OCR"
print(n_tacs)
## [1] 22
#write.csv(idps_tdcs$id, 'subjects-tdcs.csv', row.names = FALSE)
#write.csv(idps_neurocaps_mcr$id, 'subjects-neurocapsMCR.csv', row.names = FALSE)
#write.csv(idps_neurocaps_ocr$id, 'subjects-neurocapsOCR.csv', row.names = FALSE)
#write.csv(idps_tacs$id, 'subjects-tacs.csv', row.names = FALSE)


#write out behavioral data to share with Hamed/Henry

behavioral_cols <- grepl('response', names(idps_tdcs))
behavioral_cols[names(idps_tdcs) == 'id'] <- TRUE
#write.csv(idps_tdcs[, behavioral_cols], 'behavior-tdcs.csv', row.names = FALSE)

behavioral_cols <- grepl('tpre_response', names(idps_tacs))
behavioral_cols[names(idps_tacs) == 'id'] <- TRUE
#write.csv(idps_tacs[, behavioral_cols], 'behavior-tacs.csv', row.names = FALSE)

behavioral_cols <- grepl('response', names(idps_neurocaps_mcr))
behavioral_cols[names(idps_neurocaps_mcr) == 'id'] <- TRUE
#write.csv(idps_neurocaps_mcr[, behavioral_cols], 'behavior-neurocapsMCR.csv', row.names = FALSE)

behavioral_cols <- grepl('response', names(idps_neurocaps_ocr))
behavioral_cols[names(idps_neurocaps_ocr) == 'id'] <- TRUE
#write.csv(idps_neurocaps_ocr[, behavioral_cols], 'behavior-neurocapsOCR.csv', row.names = FALSE)

#idps_tacs <- read.csv('idps-tacs-abstract.csv')
#variable names are:
#for tACS
#dcr_tpre_stats_tdcsprelim_[condition].[run].0.coef_mean_[roi]
#for neurocaps
#dcr_stats_tdcsprelim_[condition].[run].0.coef_mean_[roi]

#conditions are: drug or neutral
#runs are (e.g. run 1 block 1):
#r11, r12, r13, r14

#ROIs are:
#1: VMPFC extracted at p < 0.001
#3: LSTG extracted at p < 0.001
#4: RSTG extracted at p < 0.001
#5: L Ventral Striatum extracted at p < 0.05 intersected with brainnetome ROIs 219 and 223
#6: R Ventral Striatum extracted at p < 0.05 intersected with brainnetome ROIs 220 and 224
#8: R Amygdala extracted at p < 0.05 intersected with brainnetome ROIs 212 and 214
library(reshape2)
library(lme4)
## Loading required package: Matrix
library(lmerTest)
## 
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
## 
##     lmer
## The following object is masked from 'package:stats':
## 
##     step
library(ggplot2)
library(sjstats) #for icc of mixed effects models
plot_one <- function(roi, this_label, idps, limits_spag, limits_y = c(-1, 1), prefix = 'dcr_', task = ''){
  #roi <- '1'
  #this_label <- 'VMPFC'
  
  
  #roi <- '1'
  #this_label <- 'VMPFC'
  #prefix <- 'dcr_tpre_'
  #idps <- idps_tdcs
  
  column_prefixes <- c('stats_tdcsprelim_drug.r11.0.coef_mean_',
             'stats_tdcsprelim_drug.r12.0.coef_mean_',
             'stats_tdcsprelim_drug.r13.0.coef_mean_',
             'stats_tdcsprelim_drug.r14.0.coef_mean_',
             'stats_tdcsprelim_neutral.r11.0.coef_mean_',
             'stats_tdcsprelim_neutral.r12.0.coef_mean_',
             'stats_tdcsprelim_neutral.r13.0.coef_mean_',
             'stats_tdcsprelim_neutral.r14.0.coef_mean_')
  this_roi <- c('id', 'motion', paste0(prefix, column_prefixes, roi))
  one_dataset <- idps[, this_roi]
  long_data <- melt(one_dataset, id.vars = c('id', 'motion'))
  long_data$condition <- NA
  long_data$condition[grepl('neutral', long_data$variable)] <- 'neutral'
  long_data$condition[grepl('drug', long_data$variable)] <- 'drug'
  #put neutral in the intercept
  long_data$condition <- factor(long_data$condition, levels = c('neutral', 'drug'))
  long_data$time <- NA
  long_data$time[grepl('r11', long_data$variable)] <- 1
  long_data$time[grepl('r12', long_data$variable)] <- 2
  long_data$time[grepl('r13', long_data$variable)] <- 3
  long_data$time[grepl('r14', long_data$variable)] <- 4

  #mean center on time
  long_data$time <- long_data$time - mean(long_data$time)


  #this_lme <- lmer(paste('value ~ condition * time + (1|id)'), data = long_data)
  #for checking for NA's, looks like it's all good now
  #print(long_data)
  this_lme <- lmer(paste('value ~ condition * time + (1|id/condition)'), data = long_data)
  #print(summary(this_lme))
  within_visit_iccs <- sjstats::icc(this_lme)

  n_subjects <- length(unique(long_data$id))
  means <- aggregate(long_data$value, by = list(long_data$condition, long_data$time), FUN = mean, na.rm = TRUE)
  names(means)[1:2] <- c('condition', 'time')
  sds <- aggregate(long_data$value, by = list(long_data$condition, long_data$time), FUN = sd, na.rm = TRUE)
  names(sds)[1:2] <- c('condition', 'time')


  #print(this_label)
  mean_str <- paste(this_label, '_mean', sep = '')
  names(means)[3] <- mean_str 
  

  sd_str <- paste(this_label, '_sd', sep = '')
  se_str <- paste(this_label, '_se', sep = '')
  names(sds)[3] <- sd_str
  
  plot_frame <- merge(means, sds)
  plot_frame[, se_str] <- plot_frame[, sd_str] / sqrt(n_subjects)
  plot_frame$bar_top <- plot_frame[, mean_str] + plot_frame[, se_str]
  plot_frame$bar_bottom <- plot_frame[, mean_str] - plot_frame[, se_str]
  
  plot_frame$label <- this_label
  
  ebsize <- 1.5
  linesize <- 1
  #swap condition levels back for this plot, so that drug is red and neutral is blue
  plot_frame$condition <- factor(plot_frame$condition, levels = c('drug', 'neutral'))
  p <- ggplot(plot_frame) + geom_line(aes_string(x = 'time', y = mean_str, color = 'condition'), size = linesize) + geom_errorbar(aes_string(x = 'time', ymax = 'bar_top',                                                       ymin = 'bar_bottom', color = 'condition', width = 0.1), size = ebsize) +
                      theme(text = element_text(size=10)) + ggtitle(task) + ylab('') + xlab('') + ylim(limits_y)
  
  #print(p)
  
  
  p2 <- ggplot(data = long_data[long_data$condition == 'drug',], aes_string(x = 'time', y = 'value', group = 'id')) + 
    geom_line() + stat_summary(aes(group = 1), geom = "point", fun.y = mean, shape = 17, size = 3) + 
    stat_smooth(aes(group = 1)) + labs(
    x = "Time",
    y = "% Signal Change",
    title = "drug") 
  #print(p2)
  p3 <- ggplot(data = long_data[long_data$condition == 'neutral',], aes_string(x = 'time', y = 'value', group = 'id')) + 
    geom_line() + stat_summary(aes(group = 1), geom = "point", fun.y = mean, shape = 17, size = 3) + 
    stat_smooth(aes(group = 1)) + labs(
    x = "Time",
    y = "% Signal Change",
    title = "neutral") 
  #print(p3) 
  
  wide_icc_data <- dcast(long_data, id ~ variable)
  within_iccs_simple <- ICC(wide_icc_data[, 2:ncol(wide_icc_data)])
  within_iccs_simple_intervals <- extract_icc_intervals(within_iccs_simple)
  
  
  wide_icc_data_contrasts <- dcast(long_data, id + time ~ condition)
  wide_icc_data_contrasts$DrugvNeutral <- wide_icc_data_contrasts$drug - wide_icc_data_contrasts$neutral
  contrast_icc_data <- dcast(wide_icc_data_contrasts, id ~ time, value.var = 'DrugvNeutral')
  
  contrast_iccs <- ICC(contrast_icc_data[, 2:ncol(contrast_icc_data)])
  within_contrast_icc_intervals <- extract_icc_intervals(contrast_iccs)
  
  
  return(list(model = this_lme, plotframe = plot_frame, p = p, p2=p2, p3=p3, dset = long_data, within_visit_iccs = within_visit_iccs,
              within_iccs_simple_intervals = within_iccs_simple_intervals,
              within_contrast_icc_intervals = within_contrast_icc_intervals, task = task))
  
   
}


plot_one_beh <- function(measure, this_label, idps, limits_spag, limits_y = c(1, 4), prefix = 'dcr_', task = ''){
  #same as above, but for behavioral data
  #prefix will be dcr_tpre_ for tdcs/tacs and dcr_ for neurocaps 
  #measure will be box_rt, craving_rt, or craving
  #measure <- 'craving'
  #prefix <- 'dcr_tpre_'
  #this_label <- 'Craving'
  #idps <- idps_tdcs
  
  
  #idps_tacs$dcr_tpre_response_craving_[0-7]
  #idp_tdcs
  #dcr_tpre_response_box_rt_0
  #dcr_tpre_response_craving_rt_0
  #idps_neurocaps_ocr$dcr_response_box_craving[0-7]

  column_suffixes <- 0:7
  
  these_cols <- c('id', paste0(prefix, 'response_', measure, '_', column_suffixes))
  library(reshape2)
  one_dataset <- idps[, these_cols]
  long_data <- melt(one_dataset, id.vars = c('id'))
  long_data$variable <- as.character(long_data$variable)
  long_data$condition <- NA
  long_data$number <- substr(long_data$variable, nchar(long_data$variable), nchar(long_data$variable))
  long_data$condition[long_data$number %in% c('0', '2', '4', '6')] <- 'neutral'
  long_data$condition[long_data$number %in% c('1', '3', '5', '7')] <- 'drug'
  #put neutral in the intercept
  long_data$condition <- factor(long_data$condition, levels = c('neutral', 'drug'))
  #time = block number, just like for the imaging variables
  long_data$time <- NA
  long_data$time[long_data$number %in% c('0', '1')] <- 1
  long_data$time[long_data$number %in% c('2', '3')] <- 2
  long_data$time[long_data$number %in% c('4', '5')] <- 3
  long_data$time[long_data$number %in% c('6', '7')] <- 4

  #mean center on time
  long_data$time <- long_data$time - mean(long_data$time)

  library(lme4)
  library(lmerTest)
  library(ggplot2)
  library(sjstats) #for icc of mixed effects models
  #this_lme <- lmer(paste('value ~ condition * time + (1|id)'), data = long_data)
  #for checking for NA's, looks like it's all good now
  #print(long_data)
  this_lme <- lmer(paste('value ~ condition * time +  (1|id/condition)'), data = long_data)
  #print(summary(this_lme))
  within_visit_iccs <- sjstats::icc(this_lme)

  n_subjects <- length(unique(long_data$id))
  means <- aggregate(long_data$value, by = list(long_data$condition, long_data$time), FUN = mean, na.rm = TRUE)
  names(means)[1:2] <- c('condition', 'time')
  sds <- aggregate(long_data$value, by = list(long_data$condition, long_data$time), FUN = sd, na.rm = TRUE)
  names(sds)[1:2] <- c('condition', 'time')


  #print(this_label)
  mean_str <- paste(this_label, '_mean', sep = '')
  names(means)[3] <- mean_str 
  

  sd_str <- paste(this_label, '_sd', sep = '')
  se_str <- paste(this_label, '_se', sep = '')
  names(sds)[3] <- sd_str
  
  plot_frame <- merge(means, sds)
  plot_frame[, se_str] <- plot_frame[, sd_str] / sqrt(n_subjects)
  plot_frame$bar_top <- plot_frame[, mean_str] + plot_frame[, se_str]
  plot_frame$bar_bottom <- plot_frame[, mean_str] - plot_frame[, se_str]
  
  plot_frame$label <- this_label
  
  ebsize <- 1.5
  linesize <- 1
  #swap condition levels back for this plot, so that drug is red and neutral is blue
  plot_frame$condition <- factor(plot_frame$condition, levels = c('drug', 'neutral'))
  p <- ggplot(plot_frame) + geom_line(aes_string(x = 'time', y = mean_str, color = 'condition'), size = linesize) + geom_errorbar(aes_string(x = 'time', ymax = 'bar_top',                                                       ymin = 'bar_bottom', color = 'condition', width = 0.1), size = ebsize) +
                      theme(text = element_text(size=10)) + ggtitle(task) + ylab('') + xlab('') + ylim(limits_y)
  
  #print(p)
  
  
  p2 <- ggplot(data = long_data[long_data$condition == 'drug',], aes_string(x = 'time', y = 'value', group = 'id')) + 
    geom_line() + stat_summary(aes(group = 1), geom = "point", fun.y = mean, shape = 17, size = 3) + 
    stat_smooth(aes(group = 1)) + labs(
    x = "Time",
    y = this_label,
    title = "drug") 
  #print(p2)
  p3 <- ggplot(data = long_data[long_data$condition == 'neutral',], aes_string(x = 'time', y = 'value', group = 'id')) + 
    geom_line() + stat_summary(aes(group = 1), geom = "point", fun.y = mean, shape = 17, size = 3) + 
    stat_smooth(aes(group = 1)) + labs(
    x = "Time",
    y = this_label,
    title = "neutral") 
  #print(p3) 
  
  wide_icc_data <- dcast(long_data, id ~ variable)
  within_iccs_simple <- ICC(wide_icc_data[, 2:ncol(wide_icc_data)])
  within_iccs_simple_intervals <- extract_icc_intervals(within_iccs_simple)
  
  
  wide_icc_data_contrasts <- dcast(long_data, id + time ~ condition)
  wide_icc_data_contrasts$DrugvNeutral <- wide_icc_data_contrasts$drug - wide_icc_data_contrasts$neutral
  contrast_icc_data <- dcast(wide_icc_data_contrasts, id ~ time, value.var = 'DrugvNeutral')
  
  contrast_iccs <- ICC(contrast_icc_data[, 2:ncol(contrast_icc_data)])
  within_contrast_icc_intervals <- extract_icc_intervals(contrast_iccs)
  
  
  return(list(model = this_lme, plotframe = plot_frame, p = p, p2=p2, p3=p3, dset = long_data, within_visit_iccs = within_visit_iccs,
              within_iccs_simple_intervals = within_iccs_simple_intervals,
              within_contrast_icc_intervals = within_contrast_icc_intervals, task = task))
  
   
}


library(sjPlot)
library(grid)
library(gridExtra)
library(ggpubr)
## Loading required package: magrittr
library(psych)
## 
## Attaching package: 'psych'
## The following objects are masked from 'package:sjstats':
## 
##     pca, phi
## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha
#condition by time interaction
#1: VMPFC extracted at p < 0.001
#3: LSTG extracted at p < 0.001
#4: RSTG extracted at p < 0.001
#5: L Ventral Striatum extracted at p < 0.05 intersected with brainnetome ROIs 219 and 223
#6: R Ventral Striatum extracted at p < 0.05 intersected with brainnetome ROIs 220 and 224
#8: R Amygdala extracted at p < 0.05 intersected with brainnetome ROIs 212 and 214

#main effect of condition
#9: LIFG extracted at p < 0.001
#10: RIFG at p < 0.001
#12: RDLPFC at p < 0.001

extract_icc_intervals <- function(this_icc){
  
  single_absolute <- round(this_icc$results$ICC[1], 3)
  single_absolute_lower <- this_icc$results$`lower bound`[1]
  single_absolute_upper <- this_icc$results$`upper bound`[1]
  single_fixed <- round(this_icc$results$ICC[3], 3)
  single_fixed_lower <- this_icc$results$`lower bound`[3]
  single_fixed_upper <- this_icc$results$`upper bound`[3]
  average_fixed <- round(this_icc$results$ICC[6], 3)
  average_fixed_lower <- this_icc$results$`lower bound`[6]
  average_fixed_upper <- this_icc$results$`upper bound`[6]
  this_row <- data.frame(t(c(ICC1 = single_absolute, ICC1_lower = single_absolute_lower, ICC1_upper = single_absolute_upper, ICC3 = single_fixed,
                             ICC3_lower = single_fixed_lower, ICC3_upper = single_fixed_upper,
                             ICC3k = average_fixed, ICC3k_lower = average_fixed_lower, ICC3k_upper = average_fixed_upper)))
  return(this_row)   
}

compare_models <- function(tdcs_model, neurocaps_mcr_model, neurocaps_ocr_model, tacs_model, label){
  #prints off tables of z-scores and p-values for the coefficients of condition and condition*time compared between all pairs of models
  #tdcs_summary <- summary(tdcs_model)
  #fill this with z-scores testing for differences between coefficients of different models
  ztable_condition <- data.frame(tdcs=c(NA,NA,NA,NA), neurocaps_mcr = c(NA,NA,NA,NA),
                       neurocaps_ocr = c(NA,NA,NA,NA), tacs=c(NA,NA,NA,NA))
  rownames(ztable_condition) <- c('tdcs', 'neurocaps_mcr', 'neurocaps_ocr', 'tacs')
  
  ztable_conditionbytime <- data.frame(tdcs=c(NA,NA,NA,NA), neurocaps_mcr = c(NA,NA,NA,NA),
                       neurocaps_ocr = c(NA,NA,NA,NA), tacs=c(NA,NA,NA,NA))
  rownames(ztable_conditionbytime) <- c('tdcs', 'neurocaps_mcr', 'neurocaps_ocr', 'tacs')
  
  
  
  all_models <- list(summary(tdcs_model), summary(neurocaps_mcr_model), summary(neurocaps_ocr_model), summary(tacs_model))
  for (i in 1:3){
    for (j in (i+1):4){
      #condition beta
      beta1 <- all_models[[i]]$coefficients[2,1]
      #condition SE
      se1 <- all_models[[i]]$coefficients[2,2]
      #condition beta
      beta2 <- all_models[[j]]$coefficients[2,1]
      #condition SE
      se2 <- all_models[[j]]$coefficients[2,2]
      this_z <- (beta1 - beta2) / sqrt(se1 * se1 + se2 * se2)
      ztable_condition[i,j] <- ztable_condition[j,i] <- this_z
      
      
      #conditionbytime beta--row 5 for imaging and 4 for behavioral (since motion is not included there)--but always the last row
      amodel <- all_models[[i]]$coefficients
      cbytime_idx <- nrow(amodel)
      beta1 <- all_models[[i]]$coefficients[cbytime_idx,1]
      #conditionbytime SE
      se1 <- all_models[[i]]$coefficients[cbytime_idx,2]
      #conditionbytime beta
      beta2 <- all_models[[j]]$coefficients[cbytime_idx,1]
      #conditionbytime SE
      se2 <- all_models[[j]]$coefficients[cbytime_idx,2]
      this_z <- (beta1 - beta2) / sqrt(se1 * se1 + se2 * se2)
      ztable_conditionbytime[i,j] <- ztable_conditionbytime[j,i] <- this_z
      
    }
  }
  
  ptable_condition <- ztable_condition
  ptable_conditionbytime <- ztable_conditionbytime
  for (n in names(ptable_condition)){
    ptable_condition[,n] <- 2*pnorm(-abs(ptable_condition[,n]))
    ptable_conditionbytime[,n] <- 2*pnorm(-abs(ptable_conditionbytime[,n]))
  }
  print('###Condition z-scores###')
  print(ztable_condition)
  print('###Condition by time z-scores###')
  print(ztable_conditionbytime)
    
  print('###Condition p-values###')
  print(ptable_condition)
  print('###Condition by time p-values###')
  print(ptable_conditionbytime)
  
  
  write.csv(ztable_condition, paste0('ztable-condition-', label, '.csv'), row.names = FALSE) 
  write.csv(ztable_conditionbytime, paste0('ztable-conditionbytime-', label, '.csv'), row.names = FALSE)
  
  write.csv(ptable_condition, paste0('ptable-condition-', label, '.csv'), row.names = FALSE) 
  write.csv(ptable_conditionbytime, paste0('ptable-conditionbytime-', label, '.csv'), row.names = FALSE)
  
  
  
  
}
  ##taken from jsPlot/color_utils.R to matchbe able to match colors easily
col_check2 <- function(geom.colors, collen) {
  # --------------------------------------------
  # check color argument
  # --------------------------------------------
  # check for corrct color argument
  if (!is.null(geom.colors)) {
    # check for color brewer palette
    if (is.brewer.pal(geom.colors[1])) {
      geom.colors <- scales::brewer_pal(palette = geom.colors[1])(collen)
    } else if (is.sjplot.pal(geom.colors[1])) {
      geom.colors <- get_sjplot_colorpalette(geom.colors[1], collen)
      # do we have correct amount of colours?
    } else if (geom.colors[1] == "gs") {
      geom.colors <- scales::grey_pal()(collen)
      # do we have correct amount of colours?
    } else if (geom.colors[1] == "bw") {
      geom.colors <- rep("black", times = collen)
      # do we have correct amount of colours?
    } else if (length(geom.colors) > collen) {
      # shorten palette
      geom.colors <- geom.colors[1:collen]
    } else if (length(geom.colors) < collen) {
      # repeat color palette
      geom.colors <- rep(geom.colors, times = collen)
      # shorten to required length
      geom.colors <- geom.colors[1:collen]
    }
  } else {
    geom.colors <- scales::brewer_pal(palette = "Set1")(collen)
  }

  geom.colors
}
# check whether a color value is indicating
# a color brewer palette
is.brewer.pal <- function(pal) {
  bp.seq <- c("BuGn", "BuPu", "GnBu", "OrRd", "PuBu", "PuBuGn", "PuRd", "RdPu",
              "YlGn", "YlGnBu", "YlOrBr", "YlOrRd", "Blues", "Greens", "Greys",
              "Oranges", "Purples", "Reds")
  bp.div <- c("BrBG", "PiYG", "PRGn", "PuOr", "RdBu", "RdGy", "RdYlBu",
              "RdYlGn", "Spectral")
  bp.qul <- c("Accent", "Dark2", "Paired", "Pastel1", "Pastel2", "Set1",
              "Set2", "Set3")
  bp <- c(bp.seq, bp.div, bp.qul)
  pal %in% bp
}


plot_combination <- function(roi, label, limits_spag = c(-2,2), limits_y = c(-1, 1), imaging = TRUE, forest_range = c(-0.5, 0.5)){

  #roi <- '9'
  #label <- 'LIFG'
  #limits_y = c(-0.3, 0.5)
  #limits_spag = c(-2, 2)
  #imaging = TRUE
  
  if (imaging){
    tdcs_list <- plot_one(roi, label, idps_tdcs, limits_spag = limits_spag, limits_y = limits_y, prefix = 'dcr_tpre_', task = 'Discovery Sample')
    neurocaps_mcr_list <- plot_one(roi, label, idps_neurocaps_mcr, limits_spag = limits_spag, limits_y = limits_y, prefix = 'dcr_', task = 'Replication Sample 1')
    neurocaps_ocr_list <- plot_one(roi, label, idps_neurocaps_ocr, limits_spag = limits_spag, limits_y = limits_y, prefix = 'dcr_', task = 'Replication Sample 2')
    tacs_list <- plot_one(roi, label, idps_tacs, limits_spag = limits_spag, limits_y = limits_y, prefix = 'dcr_tpre_', task = 'Sample 2-Retest')
  } else {
    tdcs_list <- plot_one_beh(roi, label, idps_tdcs, limits_spag = limits_spag, limits_y = limits_y, prefix = 'dcr_tpre_', task = 'Discovery Sample')
    neurocaps_mcr_list <- plot_one_beh(roi, label, idps_neurocaps_mcr, limits_spag = limits_spag, limits_y = limits_y, prefix = 'dcr_', task = 'Replication Sample 1')
    neurocaps_ocr_list <- plot_one_beh(roi, label, idps_neurocaps_ocr, limits_spag = limits_spag, limits_y = limits_y, prefix = 'dcr_', task = 'Replication Sample 2')
    tacs_list <- plot_one_beh(roi, label, idps_tacs, limits_spag = limits_spag, limits_y = limits_y, prefix = 'dcr_tpre_', task = 'Sample 2-Retest')
  }
  compare_models(tdcs_list$model, neurocaps_mcr_list$model, neurocaps_ocr_list$model, tacs_list$model, label)
  
  
  ###
  #extract betas and SEs for contrasts, to compare with finer grained method
  plot_frame <- NULL
  for(this_list in list(tdcs_list, neurocaps_mcr_list, neurocaps_ocr_list, tacs_list)){
    print('one')
    s1 <- summary(this_list$model)
    condition_beta <- s1$coefficients[, 'Estimate']['conditiondrug']
    conditionxtime_beta <- s1$coefficients[, 'Estimate']['conditiondrug:time']
    
    condition_se <- s1$coefficients[, 'Std. Error']['conditiondrug']
    conditionxtime_se <- s1$coefficients[, 'Std. Error']['conditiondrug:time']
    
    tr_beta <- s1$coefficients[, 'Estimate']['time']
    tr_se <- s1$coefficients[, 'Std. Error']['time']
    
    
    plot_row <- data.frame(beta = c(condition_beta, conditionxtime_beta, tr_beta), se = c(condition_se, conditionxtime_se, tr_se),
                           effect = c('drugVneutral', 'drugVneutralXtime', 'Time'), task = this_list$task)
    plot_frame <- rbind(plot_frame, plot_row)
  }
  #ggplot(plot_frame) + geom_point(aes(y = beta, x = effect))
  plot_frame$lower = plot_frame$beta - 1.96*plot_frame$se
  plot_frame$upper = plot_frame$beta + 1.96*plot_frame$se
  #plot_frame$task <- factor(plot_frame$task, levels = c('Discovery Sample', 'Replication Sample 1', 'Replication Sample 2', 'Sample 2-Retest'))
  plot_frame$task <- factor(plot_frame$task, levels = c('Sample 2-Retest', 'Replication Sample 2','Replication Sample 1', 'Discovery Sample'))
  #plot_frame$task <- factor(plot_frame$task, levels = c('1212OCR2', '1212OCR','Replication Sample 1', '1212MCR'))
  plot_frame$effect <- factor(plot_frame$effect, levels = c('drugVneutral', 'Time', 'drugVneutralXtime'))
  
  p = ggplot(data=plot_frame,
    aes(x = task,y = beta, ymin = lower, ymax = upper ))+
    geom_pointrange(aes(col=task))+
    geom_hline(aes(fill=task),yintercept =0, linetype=2)+
    xlab(label)+ ylab("Contrast Beta (95% Confidence Interval)")+
    geom_errorbar(aes(ymin=lower, ymax=upper,col=task),width=0.5,cex=1)+ 
    facet_wrap(~effect,strip.position="left",nrow=3,scales = "free_y") +
    theme(plot.title=element_text(size=16,face="bold"),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        axis.text.x=element_text(face="bold"),
        axis.title=element_text(size=12,face="bold"),
        strip.text.y = element_text(hjust=0,vjust = 1,angle=180,face="bold"))+
        coord_flip() + scale_colour_manual(values = col_check2('Set1', 4))
  print(p)
  
  #TODO correct name for this guy
  #write.csv(ptable_conditionbytime, paste0('contrasttable-', label, '.csv'), row.names = TRUE)
  ###
  
  #print model with restricted range
  print(plot_models(tdcs_list$model, neurocaps_mcr_list$model, neurocaps_ocr_list$model, tacs_list$model, m.labels = c('tDCS', 'Replication Sample 1', '1212OCR', '1212OCR2')) + ylim(forest_range))
  #plot model without restricting range so you can see where 'motion' falls
  print(plot_models(tdcs_list$model, neurocaps_mcr_list$model, neurocaps_ocr_list$model, tacs_list$model, m.labels = c('tDCS', 'Replication Sample 1', '1212OCR', '1212OCR2')))# + ylim(-0.5, 0.5))
  
  #from:
  #https://stats.stackexchange.com/questions/93540/testing-equality-of-coefficients-from-two-different-regressions
  #to test for a difference in betas, compute:
  #Z = (beta1 - beta2)/sqrt(SE(beta1)^2 + SE(beta2)^2)
  pwidth = 6
  pheight = 8
  print(ggarrange(tdcs_list$p, neurocaps_mcr_list$p, neurocaps_ocr_list$p, tacs_list$p, ncol = 2, nrow = 2, common.legend = TRUE, legend = 'bottom'))
  ggsave(paste0(label, '-lineplot.png') , ggarrange(tdcs_list$p, neurocaps_mcr_list$p, neurocaps_ocr_list$p, tacs_list$p, ncol = 2, nrow = 2, 
                                                    common.legend = TRUE, legend = 'bottom'),
         width = pwidth, height = pheight)
  print(ggarrange(tdcs_list$p2, neurocaps_mcr_list$p2, neurocaps_ocr_list$p2, tacs_list$p2, ncol = 2, nrow = 2, common.legend = TRUE, legend = 'bottom'))
  
  print(ggarrange(tdcs_list$p3, neurocaps_mcr_list$p3, neurocaps_ocr_list$p3, tacs_list$p3, ncol = 2, nrow = 2, common.legend = TRUE, legend = 'bottom'))
  
  
  #get ICCs for each of the 8 conditions
  icc_data <- merge(neurocaps_ocr_list$dset, tacs_list$dset, by = c('id', 'time', 'condition'))
  all_iccs <- NULL
  mean_values <- NA
  for (con in unique(icc_data$condition)){
    for (tim in unique(icc_data$time)){
      
      for_icc <- icc_data[icc_data$time == tim & icc_data$condition == con, c('value.x', 'value.y')]
      this_icc <- ICC(for_icc)
      this_row <- extract_icc_intervals(this_icc)
      this_row$var <- paste(tim, con)      
      all_iccs <- rbind(all_iccs, this_row)
      
      #get values for drug and neutral at this time point to put into mean_values frame
      for_means <- icc_data[icc_data$time == tim & icc_data$condition == con, c('id', 'value.x', 'value.y')]
      value_cols <- paste0(con, tim, c('NeuroCaps', 'tACS'))
      names(for_means) <- c('id', value_cols)
      mean_values <- merge(mean_values, for_means, all = TRUE)
      
    }
    #get ICCs for mean of drug and neutral
    for_icc <- icc_data[icc_data$condition == con, c('id', 'value.x', 'value.y')]
    for_icc <- aggregate(for_icc[, c('value.x', 'value.y')], by = list(for_icc$id), FUN = mean)
    value_cols <- paste0(con, c('NeuroCaps', 'tACS'))
    names(for_icc) <- c('id', value_cols)
    mean_values <- merge(mean_values, for_icc, all = TRUE)
    this_icc <- ICC(for_icc[, value_cols])
    this_row <- extract_icc_intervals(this_icc)
    this_row$var <- paste(con)      
    all_iccs <- rbind(all_iccs, this_row)
  }
  
  for (tim in unique(icc_data$time)){
    neurocaps_colname <- paste0('NeuroCapsDvN', tim)
    tacs_colname <- paste0('tACSDvN', tim)
    mean_values[, neurocaps_colname] <- mean_values[, paste0('drug', tim, 'NeuroCaps')] - mean_values[, paste0('neutral', tim, 'NeuroCaps')]
    mean_values[, tacs_colname] <- mean_values[, paste0('drug', tim, 'tACS')] - mean_values[, paste0('neutral', tim, 'tACS')]
    this_icc <- ICC(mean_values[, c(neurocaps_colname, tacs_colname)])
    this_row <- extract_icc_intervals(this_icc)
    this_row$var <- paste0(tim, 'DrugvNeutral')      
    all_iccs <- rbind(all_iccs, this_row)
  }

  #get early and late mean values
  mean_values$NeuroCapsDvNEarly <- (mean_values$`NeuroCapsDvN-1.5` + mean_values$`NeuroCapsDvN-0.5`) / 2
  mean_values$NeuroCapsDvNLate <- (mean_values$`NeuroCapsDvN1.5` + mean_values$`NeuroCapsDvN0.5`) / 2
  
  mean_values$tACSDvNEarly <- (mean_values$`tACSDvN-1.5` + mean_values$`tACSDvN-0.5`) / 2
  mean_values$tACSDvNLate <- (mean_values$`tACSDvN1.5` + mean_values$`tACSDvN0.5`) / 2
  
  #compute early and late ICCs
  this_icc <- ICC(mean_values[, c('NeuroCapsDvNEarly', 'tACSDvNEarly')])
  this_row <- extract_icc_intervals(this_icc)
  this_row$var <- paste('DrugvNeutralEarly')      
  all_iccs <- rbind(all_iccs, this_row)
  
  this_icc <- ICC(mean_values[, c('NeuroCapsDvNLate', 'tACSDvNLate')])
  this_row <- extract_icc_intervals(this_icc)
  this_row$var <- paste('DrugvNeutralLate')      
  all_iccs <- rbind(all_iccs, this_row)
  
  
  
  mean_values$NeuroCapsDvN <- mean_values$drugNeuroCaps - mean_values$neutralNeuroCaps
  mean_values$tACSDvN <- mean_values$drugtACS - mean_values$neutraltACS
  
  this_icc <- ICC(mean_values[, c('NeuroCapsDvN', 'tACSDvN')])
  this_row <- extract_icc_intervals(this_icc)
  this_row$var <- paste('DrugvNeutral')      
  all_iccs <- rbind(all_iccs, this_row)
  
  #get ICCs for mean contrast of drug and neutral in each block and in blocks 1/2 and 3/4 separately
  
  
  
  melted_iccs <- all_iccs[, c('var', 'ICC1', 'ICC1_lower', 'ICC1_upper')]
  names(melted_iccs) <- c('var', 'ICC', 'lower', 'upper')
  melted_iccs$variable <- 'ICC1'

  melted_iccs3 <- all_iccs[, c('var', 'ICC3', 'ICC3_lower', 'ICC3_upper')]
  names(melted_iccs3) <- c('var', 'ICC', 'lower', 'upper')
  melted_iccs3$variable <- 'ICC3'
  
  melted_iccs <- rbind(melted_iccs, melted_iccs3)
  #melted_iccs$var <- factor(melted_iccs$var, levels = rev(to_plot))
  melted_iccs$var <- factor(melted_iccs$var, levels = c('DrugvNeutral', 'drug', 'neutral',
                                                        "DrugvNeutralLate", "DrugvNeutralEarly",
                                                        "-1.5DrugvNeutral", "-0.5DrugvNeutral",  "0.5DrugvNeutral",   "1.5DrugvNeutral",
                                                           '1.5 drug', '0.5 drug', '-0.5 drug', '-1.5 drug',
                                                           '1.5 neutral', '0.5 neutral', '-0.5 neutral', '-1.5 neutral'))
  p <- ggplot(melted_iccs) + geom_bar(aes(x = var, fill = variable, y = ICC), stat = 'identity', position = 'dodge') + coord_flip() +
    xlab('Variable') + ylab('ICC') + ggtitle(paste0('Between Session: ', label)) + geom_errorbar(aes(x = var, fill = variable, ymin = lower, ymax = upper), stat = 'identity', position = 'dodge')
  print(p)
  print(all_iccs)
  ggsave(paste0('ICCs_', label, '.png'), p)
  
  
  #plot ICCs that come from the LME
  within_icc_results <- as.data.frame(rbind(tdcs_list$within_visit_iccs, neurocaps_mcr_list$within_visit_iccs, neurocaps_ocr_list$within_visit_iccs, tacs_list$within_visit_iccs))
  within_icc_results$dset <- c('tDCS', 'WirMCR', '1212OCR', '1212OCR2')
  within_icc_results$id_condition <- within_icc_results$id + within_icc_results$`condition:id`
  melted_within_data <- melt(within_icc_results, id.vars = c('dset', 'condition:id'))
  melted_within_data$dset <- factor(melted_within_data$dset, levels = c('1212OCR',  '1212OCR2','WirMCR', 'tDCS'))
  p <- ggplot(data = melted_within_data) + geom_bar(aes(x = dset, y = value, color = variable, fill = variable), stat = 'identity', position = 'dodge') +
    coord_flip() + ggtitle(paste0('Within Session ICCs from LME for: ', label))
  print(p)
  
  #plot within-session ICCs from psych's ICC function, treating all conditions as the same (should be similar to 'id' above)
  tdcs_list$within_iccs_simple_intervals$var <- 'tDCS'
  all_iccs <- tdcs_list$within_iccs_simple_intervals
  neurocaps_mcr_list$within_iccs_simple_intervals$var <- 'WirMCR'
  all_iccs <- rbind(all_iccs, neurocaps_mcr_list$within_iccs_simple_intervals)
  neurocaps_ocr_list$within_iccs_simple_intervals$var <- '1212OCR'
  all_iccs <- rbind(all_iccs, neurocaps_ocr_list$within_iccs_simple_intervals)
  tacs_list$within_iccs_simple_intervals$var <- '1212OCR2'
  all_iccs <- rbind(all_iccs, tacs_list$within_iccs_simple_intervals)
  
  ###plot--should functionalize this
  melted_iccs <- all_iccs[, c('var', 'ICC1', 'ICC1_lower', 'ICC1_upper')]
  names(melted_iccs) <- c('var', 'ICC', 'lower', 'upper')
  melted_iccs$variable <- 'ICC1'

  melted_iccs3 <- all_iccs[, c('var', 'ICC3', 'ICC3_lower', 'ICC3_upper')]
  names(melted_iccs3) <- c('var', 'ICC', 'lower', 'upper')
  melted_iccs3$variable <- 'ICC3'
  
  melted_iccs3k <- all_iccs[, c('var', 'ICC3k', 'ICC3k_lower', 'ICC3k_upper')]
  names(melted_iccs3k) <- c('var', 'ICC', 'lower', 'upper')
  melted_iccs3k$variable <- 'ICC3k'
  
  melted_iccs <- rbind(melted_iccs, melted_iccs3)
  melted_iccs <- rbind(melted_iccs, melted_iccs3k)
  #melted_iccs$var <- factor(melted_iccs$var, levels = rev(to_plot))
  melted_iccs$var <- factor(melted_iccs$var, levels = c('1212OCR2',  '1212OCR','WirMCR', 'tDCS'))
  p <- ggplot(melted_iccs) + geom_bar(aes(x = var, fill = variable, y = ICC), stat = 'identity', position = 'dodge') + coord_flip() +
    xlab('Variable') + ylab('ICC') + ggtitle(paste0('Within Session Simple ICC: ', label)) + geom_errorbar(aes(x = var, fill = variable, ymin = lower, ymax = upper), stat = 'identity', position = 'dodge')
  print(p)
  
  
  #plot within-session ICCs from psych's ICC function, using contrasts
  tdcs_list$within_contrast_icc_intervals$var <- 'tDCS'
  all_iccs <- tdcs_list$within_contrast_icc_intervals
  neurocaps_mcr_list$within_contrast_icc_intervals$var <- 'WirMCR'
  all_iccs <- rbind(all_iccs, neurocaps_mcr_list$within_contrast_icc_intervals)
  neurocaps_ocr_list$within_contrast_icc_intervals$var <- '1212OCR'
  all_iccs <- rbind(all_iccs, neurocaps_ocr_list$within_contrast_icc_intervals)
  tacs_list$within_contrast_icc_intervals$var <- '1212OCR2'
  all_iccs <- rbind(all_iccs, tacs_list$within_contrast_icc_intervals)
  
  ###plot--should functionalize this
  melted_iccs <- all_iccs[, c('var', 'ICC1', 'ICC1_lower', 'ICC1_upper')]
  names(melted_iccs) <- c('var', 'ICC', 'lower', 'upper')
  melted_iccs$variable <- 'ICC1'

  melted_iccs3 <- all_iccs[, c('var', 'ICC3', 'ICC3_lower', 'ICC3_upper')]
  names(melted_iccs3) <- c('var', 'ICC', 'lower', 'upper')
  melted_iccs3$variable <- 'ICC3'
  
  melted_iccs3k <- all_iccs[, c('var', 'ICC3k', 'ICC3k_lower', 'ICC3k_upper')]
  names(melted_iccs3k) <- c('var', 'ICC', 'lower', 'upper')
  melted_iccs3k$variable <- 'ICC3k'
  
  melted_iccs <- rbind(melted_iccs, melted_iccs3)
  melted_iccs <- rbind(melted_iccs, melted_iccs3k)
  #melted_iccs$var <- factor(melted_iccs$var, levels = rev(to_plot))
  melted_iccs$var <- factor(melted_iccs$var, levels = c('1212OCR2',  '1212OCR','WirMCR', 'tDCS'))
  p <- ggplot(melted_iccs) + geom_bar(aes(x = var, fill = variable, y = ICC), stat = 'identity', position = 'dodge') + coord_flip() +
    xlab('Variable') + ylab('ICC') + ggtitle(paste0('Within Session Contrast ICC: ', label)) + geom_errorbar(aes(x = var, fill = variable, ymin = lower, ymax = upper), stat = 'identity', position = 'dodge')
  print(p)
  
  ###
  
  
  
  
  print('###tDCS MCR###')
  print(summary(tdcs_list$model))
  #print(anova(tdcs_list$model, type = 'marginal'))
  print('###Neurocaps MCR###')
  print(summary(neurocaps_mcr_list$model))
  print('###Neurocaps OCR###')
  print(summary(neurocaps_ocr_list$model))
  print('###tACS OCR###')
  print(summary(tacs_list$model))
  
  
  ##make model output table to save
  tablefile <- paste0(label, '-modeltable.csv')
  write.table('tDCS Model Output', tablefile)
  write.table(round(summary(tdcs_list$model)$coefficients, digits = 3), file = tablefile,
            row.names = TRUE, append = TRUE, sep = ',')
  
  write.table('Neurocaps MCR Model Output', tablefile, append = TRUE)
  write.table(round(summary(neurocaps_mcr_list$model)$coefficients, digits = 3), file = tablefile,
            row.names = TRUE, append = TRUE, sep = ',')
  
  write.table('Neurocaps OCR Model Output', tablefile, append = TRUE)
  write.table(round(summary(neurocaps_ocr_list$model)$coefficients, digits = 3), file = tablefile,
            row.names = TRUE, append = TRUE, sep = ',')
  
  write.table('tACS OCR Model Output', tablefile, append = TRUE)
  write.table(round(summary(tacs_list$model)$coefficients, digits = 3), file = tablefile,
            row.names = TRUE, append = TRUE, sep = ',')
  
  
}

Condition

LIFG

knitr::include_graphics("LIFG_z16.png")

plot_combination('9', 'LIFG', limits_y = c(-0.15, 0.45))
## [1] "###Condition z-scores###"
##                    tdcs neurocaps_mcr neurocaps_ocr       tacs
## tdcs                 NA     1.1616086     0.8054014  0.3833179
## neurocaps_mcr 1.1616086            NA    -0.1122026 -0.5502508
## neurocaps_ocr 0.8054014    -0.1122026            NA -0.3762163
## tacs          0.3833179    -0.5502508    -0.3762163         NA
## [1] "###Condition by time z-scores###"
##                     tdcs neurocaps_mcr neurocaps_ocr      tacs
## tdcs                  NA      1.222763    -0.3421076 -1.831715
## neurocaps_mcr  1.2227631            NA    -1.2144775 -2.462810
## neurocaps_ocr -0.3421076     -1.214477            NA -1.190708
## tacs          -1.8317149     -2.462810    -1.1907079        NA
## [1] "###Condition p-values###"
##                    tdcs neurocaps_mcr neurocaps_ocr      tacs
## tdcs                 NA     0.2453945     0.4205881 0.7014841
## neurocaps_mcr 0.2453945            NA     0.9106628 0.5821474
## neurocaps_ocr 0.4205881     0.9106628            NA 0.7067561
## tacs          0.7014841     0.5821474     0.7067561        NA
## [1] "###Condition by time p-values###"
##                     tdcs neurocaps_mcr neurocaps_ocr       tacs
## tdcs                  NA    0.22141917     0.7322699 0.06699391
## neurocaps_mcr 0.22141917            NA     0.2245654 0.01378528
## neurocaps_ocr 0.73226988    0.22456544            NA 0.23376828
## tacs          0.06699391    0.01378528     0.2337683         NA
## [1] "one"
## [1] "one"
## [1] "one"
## [1] "one"
## Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

##      ICC1  ICC1_lower ICC1_upper  ICC3  ICC3_lower ICC3_upper ICC3k ICC3k_lower ICC3k_upper               var
## 1  -0.011 -0.41613414  0.4015172 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192         -0.5 drug
## 2   0.570  0.21179397  0.7944981 0.570  0.20462615  0.7956329 0.726  0.33973387   0.8861866         -1.5 drug
## 3   0.372 -0.04089754  0.6791548 0.385 -0.03340966  0.6887739 0.556 -0.06912889   0.8157088          0.5 drug
## 4   0.203 -0.22190266  0.5668081 0.203 -0.22901444  0.5688981 0.338 -0.59408230   0.7252200          1.5 drug
## 5   0.535  0.16369113  0.7753645 0.543  0.16742384  0.7810495 0.704  0.28682614   0.8770666              drug
## 6   0.060 -0.35530845  0.4597212 0.089 -0.33631783  0.4846977 0.164 -1.01349062   0.6529244      -0.5 neutral
## 7   0.383 -0.02882437  0.6856138 0.387 -0.03072059  0.6901861 0.559 -0.06338850   0.8166983      -1.5 neutral
## 8   0.069 -0.34750317  0.4667153 0.090 -0.33572472  0.4852091 0.165 -1.01079998   0.6533882       0.5 neutral
## 9   0.031 -0.38099588  0.4359664 0.066 -0.35686693  0.4666275 0.124 -1.10977631   0.6363272       1.5 neutral
## 10  0.219 -0.20686902  0.5774081 0.219 -0.21402982  0.5794607 0.359 -0.54462579   0.7337450           neutral
## 11  0.000 -0.40701472  0.4106858 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192  -0.5DrugvNeutral
## 12  0.243 -0.18220152  0.5942463 0.243 -0.18943545  0.5962381 0.391 -0.46741607   0.7470541  -1.5DrugvNeutral
## 13  0.511  0.13178270  0.7620211 0.511  0.12441301  0.7633116 0.677  0.22129414   0.8657705   0.5DrugvNeutral
## 14 -0.032 -0.43372330  0.3833580 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192   1.5DrugvNeutral
## 15  0.238 -0.18665261  0.5912574 0.238 -0.19387401  0.5932602 0.385 -0.48100177   0.7447122 DrugvNeutralEarly
## 16  0.327 -0.09261566  0.6501636 0.329 -0.09801728  0.6531131 0.495 -0.21733738   0.7901614  DrugvNeutralLate
## 17  0.491  0.10440168  0.7501278 0.491  0.09698501  0.7514737 0.658  0.17682102   0.8581045      DrugvNeutral
## Saving 7 x 5 in image

## [1] "###tDCS MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 134.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.1045 -0.5429 -0.0060  0.5391  3.4643 
## 
## Random effects:
##  Groups       Name        Variance  Std.Dev. 
##  condition:id (Intercept) 4.717e-15 6.868e-08
##  id           (Intercept) 2.123e-02 1.457e-01
##  Residual                 6.172e-02 2.484e-01
## Number of obs: 520, groups:  condition:id, 130; id, 65
## 
## Fixed effects:
##                      Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)          0.117707   0.023748 101.644762   4.956 2.88e-06 ***
## conditiondrug        0.140665   0.021789 451.999997   6.456 2.78e-10 ***
## time                 0.013784   0.013781 451.999997   1.000    0.318    
## conditiondrug:time  -0.003891   0.019489 451.999997  -0.200    0.842    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.459              
## time         0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707
## [1] "###Neurocaps MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 143.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.8143 -0.4664 -0.0174  0.5312  3.8111 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.005082 0.07129 
##  id           (Intercept) 0.021606 0.14699 
##  Residual                 0.084911 0.29139 
## Number of obs: 232, groups:  condition:id, 58; id, 29
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)  
## (Intercept)          0.04044    0.04065  46.53784   0.995   0.3250  
## conditiondrug        0.08509    0.04260  28.00000   1.998   0.0556 .
## time                 0.02114    0.02420 172.00000   0.874   0.3836  
## conditiondrug:time  -0.05205    0.03422 172.00000  -1.521   0.1301  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.524              
## time         0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707
## [1] "###Neurocaps OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 107.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7872 -0.5239  0.0764  0.5811  2.5750 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.012444 0.11155 
##  id           (Intercept) 0.006647 0.08153 
##  Residual                 0.084054 0.28992 
## Number of obs: 176, groups:  condition:id, 44; id, 22
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)  
## (Intercept)          0.10950    0.04270  40.87709   2.565   0.0141 *
## conditiondrug        0.09291    0.05515  21.00001   1.685   0.1069  
## time                 0.02024    0.02764 130.00000   0.732   0.4653  
## conditiondrug:time   0.01105    0.03909 130.00000   0.283   0.7778  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.646              
## time         0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707
## [1] "###tACS OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 120.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.1403 -0.5151 -0.0070  0.4534  2.9613 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.003373 0.05808 
##  id           (Intercept) 0.030663 0.17511 
##  Residual                 0.088214 0.29701 
## Number of obs: 176, groups:  condition:id, 44; id, 22
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)   
## (Intercept)          0.14056    0.05049  32.33604   2.784  0.00891 **
## conditiondrug        0.12043    0.04808  21.00000   2.505  0.02056 * 
## time                -0.05792    0.02832 129.99999  -2.045  0.04283 * 
## conditiondrug:time   0.07769    0.04005 129.99999   1.940  0.05455 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.476              
## time         0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707

RIFG

knitr::include_graphics("RIFG_z7.png")

plot_combination('10', 'RIFG', limits_y = c(-0.15, 0.5))
## [1] "###Condition z-scores###"
##                    tdcs neurocaps_mcr neurocaps_ocr      tacs
## tdcs                 NA     0.7123509     0.7960760 1.4727623
## neurocaps_mcr 0.7123509            NA     0.1962350 0.7743601
## neurocaps_ocr 0.7960760     0.1962350            NA 0.5133078
## tacs          1.4727623     0.7743601     0.5133078        NA
## [1] "###Condition by time z-scores###"
##                     tdcs neurocaps_mcr neurocaps_ocr      tacs
## tdcs                  NA    0.86275495    0.95744543 -1.402705
## neurocaps_mcr  0.8627549            NA    0.09483428 -1.808677
## neurocaps_ocr  0.9574454    0.09483428            NA -1.872828
## tacs          -1.4027053   -1.80867704   -1.87282850        NA
## [1] "###Condition p-values###"
##                    tdcs neurocaps_mcr neurocaps_ocr      tacs
## tdcs                 NA     0.4762475     0.4259879 0.1408151
## neurocaps_mcr 0.4762475            NA     0.8444263 0.4387179
## neurocaps_ocr 0.4259879     0.8444263            NA 0.6077361
## tacs          0.1408151     0.4387179     0.6077361        NA
## [1] "###Condition by time p-values###"
##                    tdcs neurocaps_mcr neurocaps_ocr       tacs
## tdcs                 NA    0.38827221    0.33834248 0.16070473
## neurocaps_mcr 0.3882722            NA    0.92444646 0.07050119
## neurocaps_ocr 0.3383425    0.92444646            NA 0.06109207
## tacs          0.1607047    0.07050119    0.06109207         NA
## [1] "one"
## [1] "one"
## [1] "one"
## [1] "one"
## Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

##      ICC1  ICC1_lower ICC1_upper  ICC3  ICC3_lower ICC3_upper ICC3k ICC3k_lower ICC3k_upper               var
## 1   0.133 -0.28977837  0.5154618 0.133 -0.29662688  0.5177233 0.235 -0.84344104   0.6822367         -0.5 drug
## 2   0.461  0.06586948  0.7326521 0.550  0.17634251  0.7846078 0.709  0.29981490   0.8793056         -1.5 drug
## 3   0.600  0.25478919  0.8106739 0.600  0.24776965  0.8117287 0.750  0.39714005   0.8960820          0.5 drug
## 4   0.140 -0.28298280  0.5208768 0.140 -0.28986075  0.5231210 0.246 -0.81634907   0.6869067          1.5 drug
## 5   0.472  0.08058437  0.7394312 0.472  0.07313599  0.7408262 0.641  0.13630331   0.8511202              drug
## 6   0.282 -0.14080887  0.6210504 0.282 -0.14814511  0.6229417 0.440 -0.34781773   0.7676698      -0.5 neutral
## 7   0.358 -0.05764855  0.6700045 0.358 -0.06511296  0.6717010 0.527 -0.13929590   0.8036138      -1.5 neutral
## 8   0.000 -0.40701472  0.4106858 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192       0.5 neutral
## 9  -0.045 -0.44361712  0.3728581 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192       1.5 neutral
## 10  0.063 -0.35274973  0.4620252 0.063 -0.35929266  0.4644480 0.119 -1.12155001   0.6342977           neutral
## 11  0.000 -0.40701472  0.4106858 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192  -0.5DrugvNeutral
## 12  0.249 -0.17594668  0.5984104 0.374 -0.04657415  0.6817805 0.544 -0.09769852   0.8107842  -1.5DrugvNeutral
## 13  0.608  0.26729981  0.8152265 0.608  0.26032865  0.8162586 0.756  0.41311232   0.8988352   0.5DrugvNeutral
## 14 -0.011 -0.41595961  0.4016943 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192   1.5DrugvNeutral
## 15  0.112 -0.30899658  0.4997915 0.134 -0.29592416  0.5182869 0.236 -0.84060310   0.6827259 DrugvNeutralEarly
## 16  0.650  0.33015210  0.8371177 0.650  0.32345970  0.8380384 0.788  0.48880929   0.9118835  DrugvNeutralLate
## 17  0.499  0.11567526  0.7550758 0.499  0.10827656  0.7563987 0.666  0.19539629   0.8613064      DrugvNeutral
## Saving 7 x 5 in image

## [1] "###tDCS MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 130.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2607 -0.5236  0.0328  0.6058  3.5794 
## 
## Random effects:
##  Groups       Name        Variance  Std.Dev. 
##  condition:id (Intercept) 6.670e-17 8.167e-09
##  id           (Intercept) 1.309e-02 1.144e-01
##  Residual                 6.403e-02 2.530e-01
## Number of obs: 520, groups:  condition:id, 130; id, 65
## 
## Fixed effects:
##                      Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)         6.582e-02  2.116e-02  1.193e+02   3.111  0.00234 ** 
## conditiondrug       1.661e-01  2.219e-02  4.520e+02   7.485 3.76e-13 ***
## time                4.917e-03  1.404e-02  4.520e+02   0.350  0.72625    
## conditiondrug:time -6.029e-04  1.985e-02  4.520e+02  -0.030  0.97578    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.524              
## time         0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707
## [1] "###Neurocaps MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 180.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.0410 -0.3870 -0.0026  0.4170  5.8072 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.00122  0.03493 
##  id           (Intercept) 0.02198  0.14826 
##  Residual                 0.10432  0.32299 
## Number of obs: 232, groups:  condition:id, 58; id, 29
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)          0.14614    0.04122  46.70774   3.545 0.000904 ***
## conditiondrug        0.13141    0.04339  28.00000   3.029 0.005233 ** 
## time                 0.01599    0.02682 172.00001   0.596 0.551939    
## conditiondrug:time  -0.03754    0.03793 172.00001  -0.990 0.323738    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.526              
## time         0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707
## [1] "###Neurocaps OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 107.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4704 -0.5184  0.0064  0.5092  4.3247 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.014729 0.12136 
##  id           (Intercept) 0.003274 0.05722 
##  Residual                 0.084736 0.29110 
## Number of obs: 176, groups:  condition:id, 44; id, 22
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)  
## (Intercept)          0.10868    0.04220  41.70892   2.575   0.0137 *
## conditiondrug        0.11733    0.05714  21.00001   2.053   0.0527 .
## time                 0.02199    0.02775 130.00000   0.792   0.4297  
## conditiondrug:time  -0.04272    0.03925 130.00000  -1.088   0.2785  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.677              
## time         0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707
## [1] "###tACS OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 114.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5570 -0.4614  0.0175  0.4512  4.3707 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.01733  0.1317  
##  id           (Intercept) 0.04960  0.2227  
##  Residual                 0.07431  0.2726  
## Number of obs: 176, groups:  condition:id, 44; id, 22
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)
## (Intercept)          0.08769    0.06234  31.42680   1.407    0.169
## conditiondrug        0.07586    0.05714  21.00000   1.328    0.199
## time                -0.03981    0.02599 130.00000  -1.532    0.128
## conditiondrug:time   0.05800    0.03676 130.00000   1.578    0.117
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.458              
## time         0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707

RDLPFC

knitr::include_graphics("RDLPFC_z7.png")

plot_combination('12', 'RDLPFC', limits_y = c(-0.45, 0.1))
## [1] "###Condition z-scores###"
##                      tdcs neurocaps_mcr neurocaps_ocr        tacs
## tdcs                   NA    0.31425130   0.009064596  0.21861557
## neurocaps_mcr 0.314251300            NA  -0.212396545 -0.05095833
## neurocaps_ocr 0.009064596   -0.21239655            NA  0.15452835
## tacs          0.218615569   -0.05095833   0.154528347          NA
## [1] "###Condition by time z-scores###"
##                      tdcs neurocaps_mcr neurocaps_ocr        tacs
## tdcs                   NA    0.05074138    0.02546767 -0.07062211
## neurocaps_mcr  0.05074138            NA   -0.02233374 -0.09538805
## neurocaps_ocr  0.02546767   -0.02233374            NA -0.07761132
## tacs          -0.07062211   -0.09538805   -0.07761132          NA
## [1] "###Condition p-values###"
##                    tdcs neurocaps_mcr neurocaps_ocr      tacs
## tdcs                 NA     0.7533302     0.9927676 0.8269495
## neurocaps_mcr 0.7533302            NA     0.8317977 0.9593587
## neurocaps_ocr 0.9927676     0.8317977            NA 0.8771932
## tacs          0.8269495     0.9593587     0.8771932        NA
## [1] "###Condition by time p-values###"
##                    tdcs neurocaps_mcr neurocaps_ocr      tacs
## tdcs                 NA     0.9595316     0.9796819 0.9436985
## neurocaps_mcr 0.9595316            NA     0.9821817 0.9240066
## neurocaps_ocr 0.9796819     0.9821817            NA 0.9381372
## tacs          0.9436985     0.9240066     0.9381372        NA
## [1] "one"
## [1] "one"
## [1] "one"
## [1] "one"
## Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

##     ICC1  ICC1_lower ICC1_upper  ICC3  ICC3_lower ICC3_upper ICC3k ICC3k_lower ICC3k_upper               var
## 1  0.037 -0.37551627  0.4411311 0.037 -0.38193421  0.4436122 0.072 -1.23590147   0.6145864         -0.5 drug
## 2  0.153 -0.27069876  0.5305039 0.202 -0.23015741  0.5680815 0.337 -0.59793370   0.7245561         -1.5 drug
## 3  0.071 -0.34593490  0.4681084 0.071 -0.35251375  0.4705138 0.133 -1.08886868   0.6399311          0.5 drug
## 4  0.287 -0.13637989  0.6238166 0.300 -0.12966353  0.6343391 0.461 -0.29796184   0.7762638          1.5 drug
## 5  0.236 -0.18892048  0.5897263 0.236 -0.19613539  0.5917346 0.382 -0.48798116   0.7435092              drug
## 6  0.000 -0.40701472  0.4106858 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192      -0.5 neutral
## 7  0.312 -0.10857039  0.6407610 0.357 -0.06553960  0.6714658 0.527 -0.14027262   0.8034454      -1.5 neutral
## 8  0.000 -0.40701472  0.4106858 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192       0.5 neutral
## 9  0.000 -0.40701472  0.4106858 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192       1.5 neutral
## 10 0.000 -0.40701472  0.4106858 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192           neutral
## 11 0.000 -0.40701472  0.4106858 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192  -0.5DrugvNeutral
## 12 0.508  0.12668409  0.7598382 0.508  0.11930479  0.7611390 0.673  0.21317660   0.8643713  -1.5DrugvNeutral
## 13 0.177 -0.24786607  0.5478653 0.177 -0.25488525  0.5500205 0.301 -0.68415032   0.7096945   0.5DrugvNeutral
## 14 0.000 -0.40701472  0.4106858 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192   1.5DrugvNeutral
## 15 0.044 -0.36929027  0.4469343 0.044 -0.37574316  0.4493995 0.085 -1.20380949   0.6201182 DrugvNeutralEarly
## 16 0.385 -0.02560338  0.6873182 0.385 -0.03308957  0.6889423 0.556 -0.06844391   0.8158269  DrugvNeutralLate
## 17 0.305 -0.11660435  0.6359398 0.305 -0.12398856  0.6377734 0.467 -0.28307521   0.7788298      DrugvNeutral
## Saving 7 x 5 in image

## [1] "###tDCS MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 270
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.4695 -0.5557  0.0332  0.5570  4.0540 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.002227 0.04719 
##  id           (Intercept) 0.019266 0.13880 
##  Residual                 0.081454 0.28540 
## Number of obs: 520, groups:  condition:id, 130; id, 65
## 
## Fixed effects:
##                      Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)         -0.063033   0.025376 105.622506  -2.484   0.0146 *  
## conditiondrug       -0.147297   0.026365  64.000071  -5.587 5.08e-07 ***
## time                -0.037129   0.015831 388.000112  -2.345   0.0195 *  
## conditiondrug:time  -0.009138   0.022389 388.000112  -0.408   0.6834    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.519              
## time         0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707
## [1] "###Neurocaps MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 266.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.3603 -0.3804  0.0355  0.3806  4.7200 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.00000  0.0000  
##  id           (Intercept) 0.03132  0.1770  
##  Residual                 0.15362  0.3919  
## Number of obs: 232, groups:  condition:id, 58; id, 29
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)   
## (Intercept)         -0.07738    0.04903  52.27281  -1.578  0.12057   
## conditiondrug       -0.16547    0.05146 200.00000  -3.215  0.00152 **
## time                -0.01313    0.03255 200.00000  -0.403  0.68713   
## conditiondrug:time  -0.01174    0.04603 200.00000  -0.255  0.79902   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.525              
## time         0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707
## [1] "###Neurocaps OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 141.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8505 -0.4083  0.0700  0.3988  4.0875 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.02182  0.1477  
##  id           (Intercept) 0.01552  0.1246  
##  Residual                 0.09622  0.3102  
## Number of obs: 176, groups:  condition:id, 44; id, 22
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)  
## (Intercept)         -0.02715    0.05282  39.47828  -0.514   0.6102  
## conditiondrug       -0.14793    0.06458  21.00000  -2.291   0.0324 *
## time                -0.02566    0.02958 130.00000  -0.868   0.3871  
## conditiondrug:time  -0.01035    0.04183 130.00000  -0.247   0.8050  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.611              
## time         0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707
## [1] "###tACS OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 169.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.9009 -0.3164  0.0318  0.4854  3.4737 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.008613 0.09281 
##  id           (Intercept) 0.023749 0.15411 
##  Residual                 0.119896 0.34626 
## Number of obs: 176, groups:  condition:id, 44; id, 22
## 
## Fixed effects:
##                      Estimate Std. Error         df t value Pr(>|t|)  
## (Intercept)         -0.082844   0.053230  36.676582  -1.556   0.1282  
## conditiondrug       -0.161470   0.059228  20.999999  -2.726   0.0127 *
## time                -0.016560   0.033015 129.999996  -0.502   0.6168  
## conditiondrug:time  -0.005481   0.046690 129.999996  -0.117   0.9067  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.556              
## time         0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707

ConditionXTime

VMPFC

knitr::include_graphics("VMPFC_zm11.png")

plot_combination('1', 'VMPFC', limits_y = c(-0.4, 0.5))
## [1] "###Condition z-scores###"
##                     tdcs neurocaps_mcr neurocaps_ocr     tacs
## tdcs                  NA    0.07598867     0.5455601 2.524548
## neurocaps_mcr 0.07598867            NA     0.4319498 2.190228
## neurocaps_ocr 0.54556014    0.43194983            NA 1.565483
## tacs          2.52454834    2.19022821     1.5654825       NA
## [1] "###Condition by time z-scores###"
##                     tdcs neurocaps_mcr neurocaps_ocr      tacs
## tdcs                  NA    -0.7157979     0.4934780 -2.132664
## neurocaps_mcr -0.7157979            NA     0.9771868 -1.340874
## neurocaps_ocr  0.4934780     0.9771868            NA -2.096726
## tacs          -2.1326637    -1.3408743    -2.0967259        NA
## [1] "###Condition p-values###"
##                     tdcs neurocaps_mcr neurocaps_ocr       tacs
## tdcs                  NA    0.93942811     0.5853683 0.01158471
## neurocaps_mcr 0.93942811            NA     0.6657779 0.02850769
## neurocaps_ocr 0.58536833    0.66577788            NA 0.11746981
## tacs          0.01158471    0.02850769     0.1174698         NA
## [1] "###Condition by time p-values###"
##                     tdcs neurocaps_mcr neurocaps_ocr       tacs
## tdcs                  NA     0.4741162    0.62167491 0.03295232
## neurocaps_mcr 0.47411616            NA    0.32847668 0.17996127
## neurocaps_ocr 0.62167491     0.3284767            NA 0.03601785
## tacs          0.03295232     0.1799613    0.03601785         NA
## [1] "one"
## [1] "one"
## [1] "one"
## [1] "one"
## Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

##      ICC1  ICC1_lower ICC1_upper  ICC3  ICC3_lower ICC3_upper ICC3k ICC3k_lower ICC3k_upper               var
## 1   0.073 -0.34407149  0.4697585 0.073 -0.35066003  0.4721591 0.136 -1.08005068   0.6414511         -0.5 drug
## 2   0.079 -0.33889315  0.4743142 0.223 -0.20974246  0.5824346 0.364 -0.53082051   0.7361247         -1.5 drug
## 3   0.042 -0.37171705  0.4446805 0.042 -0.37815638  0.4471519 0.080 -1.21624271   0.6179750          0.5 drug
## 4   0.000 -0.40701472  0.4106858 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192          1.5 drug
## 5   0.114 -0.30707170  0.5013852 0.114 -0.31384216  0.5036910 0.205 -0.91478125   0.6699395              drug
## 6  -0.021 -0.42418129  0.3932883 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192      -0.5 neutral
## 7   0.334 -0.08416782  0.6550517 0.345 -0.08013093  0.6633308 0.512 -0.17422246   0.7975933      -1.5 neutral
## 8  -0.056 -0.45308572  0.3626098 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192       0.5 neutral
## 9   0.000 -0.40701472  0.4106858 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192       1.5 neutral
## 10 -0.037 -0.43721342  0.3796780 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192           neutral
## 11 -0.001 -0.40777027  0.4099324 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192  -0.5DrugvNeutral
## 12  0.121 -0.30082855  0.5065169 0.314 -0.11408880  0.6436927 0.478 -0.25756262   0.7832276  -1.5DrugvNeutral
## 13 -0.002 -0.40857976  0.4091240 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192   0.5DrugvNeutral
## 14  0.431  0.02938104  0.7152510 0.448  0.04316295  0.7269512 0.619  0.08275399   0.8418897   1.5DrugvNeutral
## 15 -0.048 -0.44653430  0.3697218 0.070 -0.35369437  0.4694629 0.130 -1.09451116   0.6389585 DrugvNeutralEarly
## 16  0.227 -0.19843003  0.5832447 0.227 -0.20561686  0.5852764 0.370 -0.51767681   0.7383904  DrugvNeutralLate
## 17  0.234 -0.19073916  0.5884945 0.311 -0.11713950  0.6418783 0.475 -0.26536355   0.7818829      DrugvNeutral
## Saving 7 x 5 in image

## [1] "###tDCS MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 764
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.4353 -0.4896 -0.0141  0.5046  5.2473 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.003389 0.05821 
##  id           (Intercept) 0.040859 0.20214 
##  Residual                 0.217516 0.46639 
## Number of obs: 520, groups:  condition:id, 130; id, 65
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)         -0.04460    0.03895 109.24965  -1.145   0.2548    
## conditiondrug        0.19212    0.04216  64.00003   4.557 2.39e-05 ***
## time                 0.05140    0.02587 388.00004   1.987   0.0477 *  
## conditiondrug:time  -0.18186    0.03659 388.00004  -4.971 1.00e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.541              
## time         0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707
## [1] "###Neurocaps MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 317.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1094 -0.5469  0.0178  0.5221  3.4627 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.000841 0.0290  
##  id           (Intercept) 0.012000 0.1095  
##  Residual                 0.205415 0.4532  
## Number of obs: 232, groups:  condition:id, 58; id, 29
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)   
## (Intercept)         -0.05185    0.04705  54.10927  -1.102  0.27530   
## conditiondrug        0.18654    0.06000  28.00013   3.109  0.00428 **
## time                 0.08372    0.03764 172.00014   2.224  0.02743 * 
## conditiondrug:time  -0.13563    0.05323 172.00014  -2.548  0.01171 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.638              
## time         0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707
## [1] "###Neurocaps OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 282.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0596 -0.5953  0.0095  0.5922  3.7363 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.00000  0.0000  
##  id           (Intercept) 0.02454  0.1567  
##  Residual                 0.25327  0.5033  
## Number of obs: 176, groups:  condition:id, 44; id, 22
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)   
## (Intercept)         -0.08044    0.06319  49.15502  -1.273  0.20902   
## conditiondrug        0.14476    0.07587 151.00000   1.908  0.05828 . 
## time                 0.06562    0.04798 151.00000   1.367  0.17351   
## conditiondrug:time  -0.21991    0.06786 151.00000  -3.241  0.00147 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.600              
## time         0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707
## [1] "###tACS OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 269.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1579 -0.5297  0.0025  0.5317  4.3813 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.00000  0.0000  
##  id           (Intercept) 0.02972  0.1724  
##  Residual                 0.23074  0.4804  
## Number of obs: 176, groups:  condition:id, 44; id, 22
## 
## Fixed effects:
##                      Estimate Std. Error         df t value Pr(>|t|)
## (Intercept)          0.053982   0.063033  45.250379   0.856    0.396
## conditiondrug       -0.019427   0.072416 151.000000  -0.268    0.789
## time                 0.006036   0.045800 151.000000   0.132    0.895
## conditiondrug:time  -0.023216   0.064770 151.000000  -0.358    0.721
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.574              
## time         0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707

LSTG

knitr::include_graphics("LSFG_z11.png")

plot_combination('3', 'LSTG', limits_y = c(-0.5, 0.3))
## [1] "###Condition z-scores###"
##                   tdcs neurocaps_mcr neurocaps_ocr       tacs
## tdcs                NA     1.9347464     1.0097603  1.2324662
## neurocaps_mcr 1.934746            NA    -0.7885670 -0.6726636
## neurocaps_ocr 1.009760    -0.7885670            NA  0.1405627
## tacs          1.232466    -0.6726636     0.1405627         NA
## [1] "###Condition by time z-scores###"
##                     tdcs neurocaps_mcr neurocaps_ocr      tacs
## tdcs                  NA    -0.7002844     0.3893853 -2.225859
## neurocaps_mcr -0.7002844            NA     0.8720940 -1.115811
## neurocaps_ocr  0.3893853     0.8720940            NA -2.073894
## tacs          -2.2258591    -1.1158114    -2.0738943        NA
## [1] "###Condition p-values###"
##                     tdcs neurocaps_mcr neurocaps_ocr      tacs
## tdcs                  NA    0.05302142     0.3126101 0.2177750
## neurocaps_mcr 0.05302142            NA     0.4303651 0.5011613
## neurocaps_ocr 0.31261014    0.43036512            NA 0.8882154
## tacs          0.21777498    0.50116133     0.8882154        NA
## [1] "###Condition by time p-values###"
##                     tdcs neurocaps_mcr neurocaps_ocr       tacs
## tdcs                  NA     0.4837497    0.69699112 0.02602362
## neurocaps_mcr 0.48374969            NA    0.38315709 0.26450288
## neurocaps_ocr 0.69699112     0.3831571            NA 0.03808913
## tacs          0.02602362     0.2645029    0.03808913         NA
## [1] "one"
## [1] "one"
## [1] "one"
## [1] "one"
## Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

##      ICC1 ICC1_lower ICC1_upper  ICC3 ICC3_lower ICC3_upper ICC3k ICC3k_lower ICC3k_upper               var
## 1   0.000 -0.4070147  0.4106858 0.000 -0.4132470  0.4132470 0.000  -1.4085895   0.5848192         -0.5 drug
## 2   0.122 -0.3002244  0.5070105 0.146 -0.2846024  0.5272714 0.254  -0.7956483   0.6904750         -1.5 drug
## 3   0.215 -0.2106148  0.5747915 0.219 -0.2138283  0.5796010 0.359  -0.5439736   0.7338575          0.5 drug
## 4   0.000 -0.4070147  0.4106858 0.000 -0.4132470  0.4132470 0.000  -1.4085895   0.5848192          1.5 drug
## 5   0.078 -0.3394387  0.4738363 0.078 -0.3460512  0.4762250 0.145  -1.0583433   0.6451929              drug
## 6   0.000 -0.4070147  0.4106858 0.000 -0.4132470  0.4132470 0.000  -1.4085895   0.5848192      -0.5 neutral
## 7   0.135 -0.2874377  0.5173342 0.217 -0.2156613  0.5783235 0.357  -0.5499188   0.7328326      -1.5 neutral
## 8  -0.002 -0.4084223  0.4092813 0.000 -0.4132470  0.4132470 0.000  -1.4085895   0.5848192       0.5 neutral
## 9   0.038 -0.3747622  0.4418376 0.055 -0.3664769  0.4579340 0.104  -1.1569490   0.6281958       1.5 neutral
## 10  0.000 -0.4070147  0.4106858 0.000 -0.4132470  0.4132470 0.000  -1.4085895   0.5848192           neutral
## 11  0.000 -0.4070147  0.4106858 0.000 -0.4132470  0.4132470 0.000  -1.4085895   0.5848192  -0.5DrugvNeutral
## 12 -0.009 -0.4148260  0.4028425 0.126 -0.3030504  0.5125382 0.224  -0.8696480   0.6777193  -1.5DrugvNeutral
## 13  0.113 -0.3082963  0.5003719 0.113 -0.3150611  0.5026808 0.203  -0.9199683   0.6690453   0.5DrugvNeutral
## 14 -0.003 -0.4096944  0.4080088 0.007 -0.4070884  0.4193681 0.015  -1.3731840   0.5909223   1.5DrugvNeutral
## 15 -0.009 -0.4142570  0.4034180 0.000 -0.4132470  0.4132470 0.000  -1.4085895   0.5848192 DrugvNeutralEarly
## 16  0.110 -0.3107377  0.4983452 0.110 -0.3174910  0.5006603 0.199  -0.9303644   0.6672533  DrugvNeutralLate
## 17  0.282 -0.1409204  0.6209804 0.282 -0.1482564  0.6228720 0.440  -0.3481246   0.7676169      DrugvNeutral
## Saving 7 x 5 in image

## [1] "###tDCS MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 491.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.6714 -0.5129  0.0291  0.5248  4.7424 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.00000  0.0000  
##  id           (Intercept) 0.01809  0.1345  
##  Residual                 0.13263  0.3642  
## Number of obs: 520, groups:  condition:id, 130; id, 65
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)         -0.09040    0.02808 135.46160  -3.220  0.00161 ** 
## conditiondrug       -0.04399    0.03194 452.00000  -1.377  0.16914    
## time                 0.10007    0.02020 452.00000   4.954 1.03e-06 ***
## conditiondrug:time  -0.15242    0.02857 452.00000  -5.335 1.51e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.569              
## time         0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707
## [1] "###Neurocaps MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 358.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.3322 -0.4273 -0.0179  0.3480  8.7958 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.00000  0.0000  
##  id           (Intercept) 0.03375  0.1837  
##  Residual                 0.23581  0.4856  
## Number of obs: 232, groups:  condition:id, 58; id, 29
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)   
## (Intercept)          0.01396    0.05654  58.41679   0.247  0.80590   
## conditiondrug       -0.18197    0.06376 200.00000  -2.854  0.00477 **
## time                 0.07866    0.04033 200.00000   1.951  0.05250 . 
## conditiondrug:time  -0.10776    0.05703 200.00000  -1.889  0.06028 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.564              
## time         0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707
## [1] "###Neurocaps OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 205.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7491 -0.4245  0.0577  0.4919  3.8096 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.00000  0.0000  
##  id           (Intercept) 0.02027  0.1424  
##  Residual                 0.15951  0.3994  
## Number of obs: 176, groups:  condition:id, 44; id, 22
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)         -0.02341    0.05229  45.43395  -0.448 0.656469    
## conditiondrug       -0.11281    0.06021 151.00000  -1.874 0.062913 .  
## time                 0.13143    0.03808 151.00000   3.451 0.000723 ***
## conditiondrug:time  -0.17616    0.05385 151.00000  -3.271 0.001327 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.576              
## time         0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707
## [1] "###tACS OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 204.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.2469 -0.4307  0.0228  0.4452  3.5914 
## 
## Random effects:
##  Groups       Name        Variance  Std.Dev. 
##  condition:id (Intercept) 7.375e-18 2.716e-09
##  id           (Intercept) 6.649e-02 2.579e-01
##  Residual                 1.427e-01 3.777e-01
## Number of obs: 176, groups:  condition:id, 44; id, 22
## 
## Fixed effects:
##                      Estimate Std. Error         df t value Pr(>|t|)  
## (Intercept)         -0.008088   0.068146  30.632068  -0.119   0.9063  
## conditiondrug       -0.124459   0.056947 151.000000  -2.186   0.0304 *
## time                 0.014761   0.036016 151.000000   0.410   0.6825  
## conditiondrug:time  -0.022435   0.050935 151.000000  -0.440   0.6602  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.418              
## time         0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707

RSTG

knitr::include_graphics("RSFG_z7.png")

plot_combination('4', 'RSTG', limits_y = c(-0.4, 0.25))
## [1] "###Condition z-scores###"
##                    tdcs neurocaps_mcr neurocaps_ocr       tacs
## tdcs                 NA     2.2100764     0.3935193  1.6060658
## neurocaps_mcr 2.2100764            NA    -1.5704736 -0.2390023
## neurocaps_ocr 0.3935193    -1.5704736            NA  1.1447217
## tacs          1.6060658    -0.2390023     1.1447217         NA
## [1] "###Condition by time z-scores###"
##                     tdcs neurocaps_mcr neurocaps_ocr       tacs
## tdcs                  NA    -1.3968479    -0.6932403 -2.4421241
## neurocaps_mcr -1.3968479            NA     0.5854781 -0.9980761
## neurocaps_ocr -0.6932403     0.5854781            NA -1.5448472
## tacs          -2.4421241    -0.9980761    -1.5448472         NA
## [1] "###Condition p-values###"
##                     tdcs neurocaps_mcr neurocaps_ocr      tacs
## tdcs                  NA    0.02709986     0.6939360 0.1082595
## neurocaps_mcr 0.02709986            NA     0.1163050 0.8111038
## neurocaps_ocr 0.69393599    0.11630498            NA 0.2523245
## tacs          0.10825946    0.81110379     0.2523245        NA
## [1] "###Condition by time p-values###"
##                     tdcs neurocaps_mcr neurocaps_ocr       tacs
## tdcs                  NA     0.1624593     0.4881588 0.01460112
## neurocaps_mcr 0.16245931            NA     0.5582263 0.31824246
## neurocaps_ocr 0.48815879     0.5582263            NA 0.12238323
## tacs          0.01460112     0.3182425     0.1223832         NA
## [1] "one"
## [1] "one"
## [1] "one"
## [1] "one"
## Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

##      ICC1  ICC1_lower ICC1_upper  ICC3  ICC3_lower ICC3_upper ICC3k ICC3k_lower ICC3k_upper               var
## 1   0.000 -0.40701472  0.4106858 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192         -0.5 drug
## 2   0.395 -0.01412295  0.6933300 0.395 -0.02161320  0.6949285 0.566 -0.04418130   0.8200092         -1.5 drug
## 3   0.386 -0.02518190  0.6875407 0.386 -0.03266827  0.6891638 0.557 -0.06754306   0.8159822          0.5 drug
## 4   0.000 -0.40701472  0.4106858 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192          1.5 drug
## 5   0.165 -0.25966198  0.5389808 0.165 -0.26663577  0.5411658 0.283 -0.72715781   0.7022811              drug
## 6   0.067 -0.34948960  0.4649449 0.067 -0.35604980  0.4673594 0.126 -1.10583024   0.6370074      -0.5 neutral
## 7   0.070 -0.34684083  0.4673042 0.090 -0.33541924  0.4854722 0.166 -1.00941603   0.6536268      -1.5 neutral
## 8   0.000 -0.40701472  0.4106858 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192       0.5 neutral
## 9   0.034 -0.37789252  0.4388980 0.036 -0.38257773  0.4430066 0.070 -1.23927418   0.6140050       1.5 neutral
## 10  0.000 -0.40701472  0.4106858 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192           neutral
## 11 -0.012 -0.41688235  0.4007577 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192  -0.5DrugvNeutral
## 12  0.116 -0.30546157  0.5027142 0.166 -0.26578608  0.5418121 0.284 -0.72400175   0.7028251  -1.5DrugvNeutral
## 13  0.193 -0.23183343  0.5596598 0.193 -0.23891098  0.5617747 0.324 -0.62781349   0.7194056   0.5DrugvNeutral
## 14  0.050 -0.36450590  0.4513477 0.050 -0.37098524  0.4538006 0.095 -1.17957565   0.6242955   1.5DrugvNeutral
## 15 -0.003 -0.40969686  0.4080063 0.050 -0.37129867  0.4535119 0.094 -1.18116075   0.6240223 DrugvNeutralEarly
## 16  0.167 -0.25744436  0.5406647 0.167 -0.26442684  0.5428441 0.286 -0.71896816   0.7036928  DrugvNeutralLate
## 17  0.283 -0.13993787  0.6215959 0.304 -0.12516779  0.6370622 0.466 -0.28615267   0.7782994      DrugvNeutral
## Saving 7 x 5 in image

## [1] "###tDCS MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 462
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.2283 -0.4786  0.0226  0.5172  4.4717 
## 
## Random effects:
##  Groups       Name        Variance  Std.Dev. 
##  condition:id (Intercept) 3.728e-15 6.105e-08
##  id           (Intercept) 2.009e-02 1.417e-01
##  Residual                 1.237e-01 3.517e-01
## Number of obs: 520, groups:  condition:id, 130; id, 65
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)         -0.04823    0.02802 128.34071  -1.721   0.0876 .  
## conditiondrug       -0.03724    0.03085 452.00000  -1.207   0.2281    
## time                 0.09463    0.01951 452.00000   4.850 1.70e-06 ***
## conditiondrug:time  -0.15258    0.02759 452.00000  -5.530 5.43e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.551              
## time         0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707
## [1] "###Neurocaps MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 242.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.4950 -0.3987  0.0031  0.4738  5.6161 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.006746 0.08213 
##  id           (Intercept) 0.030134 0.17359 
##  Residual                 0.132720 0.36431 
## Number of obs: 232, groups:  condition:id, 58; id, 29
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)   
## (Intercept)         -0.02264    0.04915  47.25719  -0.461  0.64713   
## conditiondrug       -0.17176    0.05247  28.00000  -3.273  0.00283 **
## time                 0.05542    0.03025 172.00000   1.832  0.06869 . 
## conditiondrug:time  -0.08147    0.04279 172.00000  -1.904  0.05857 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.534              
## time         0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707
## [1] "###Neurocaps OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 131.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8655 -0.5059  0.0564  0.5771  4.1954 
## 
## Random effects:
##  Groups       Name        Variance  Std.Dev. 
##  condition:id (Intercept) 8.943e-16 2.990e-08
##  id           (Intercept) 1.572e-02 1.254e-01
##  Residual                 1.026e-01 3.202e-01
## Number of obs: 176, groups:  condition:id, 44; id, 22
## 
## Fixed effects:
##                      Estimate Std. Error         df t value Pr(>|t|)   
## (Intercept)         -0.005619   0.043359  42.899148  -0.130  0.89749   
## conditiondrug       -0.059782   0.048279 151.000001  -1.238  0.21754   
## time                 0.091057   0.030534 151.000001   2.982  0.00334 **
## conditiondrug:time  -0.117059   0.043182 151.000001  -2.711  0.00749 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.557              
## time         0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707
## [1] "###tACS OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 190
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.1373 -0.4562  0.0596  0.4561  4.2749 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.01385  0.1177  
##  id           (Intercept) 0.04615  0.2148  
##  Residual                 0.12704  0.3564  
## Number of obs: 176, groups:  condition:id, 44; id, 22
## 
## Fixed effects:
##                      Estimate Std. Error         df t value Pr(>|t|)  
## (Intercept)          0.030926   0.064584  33.520398   0.479   0.6352  
## conditiondrug       -0.151911   0.064393  21.000000  -2.359   0.0281 *
## time                -0.007535   0.033984 129.999999  -0.222   0.8249  
## conditiondrug:time  -0.017246   0.048061 129.999999  -0.359   0.7203  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.499              
## time         0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707

LVStriatum

knitr::include_graphics("LVStri_zm1.png")

plot_combination('5', 'LVSTRI', limits_y = c(-0.22, 0.35))
## [1] "###Condition z-scores###"
##                    tdcs neurocaps_mcr neurocaps_ocr     tacs
## tdcs                 NA     0.7779924     0.1985035 2.713589
## neurocaps_mcr 0.7779924            NA    -0.3694164 1.641814
## neurocaps_ocr 0.1985035    -0.3694164            NA 1.762054
## tacs          2.7135890     1.6418139     1.7620536       NA
## [1] "###Condition by time z-scores###"
##                     tdcs neurocaps_mcr neurocaps_ocr      tacs
## tdcs                  NA     0.1731414     0.9808995 -2.384295
## neurocaps_mcr  0.1731414            NA     0.7411386 -2.132148
## neurocaps_ocr  0.9808995     0.7411386            NA -2.544832
## tacs          -2.3842952    -2.1321485    -2.5448317        NA
## [1] "###Condition p-values###"
##                      tdcs neurocaps_mcr neurocaps_ocr        tacs
## tdcs                   NA     0.4365735    0.84265117 0.006655868
## neurocaps_mcr 0.436573517            NA    0.71181736 0.100628575
## neurocaps_ocr 0.842651170     0.7118174            NA 0.078060239
## tacs          0.006655868     0.1006286    0.07806024          NA
## [1] "###Condition by time p-values###"
##                     tdcs neurocaps_mcr neurocaps_ocr       tacs
## tdcs                  NA    0.86254031    0.32664230 0.01711187
## neurocaps_mcr 0.86254031            NA    0.45860940 0.03299464
## neurocaps_ocr 0.32664230    0.45860940            NA 0.01093304
## tacs          0.01711187    0.03299464    0.01093304         NA
## [1] "one"
## [1] "one"
## [1] "one"
## [1] "one"
## Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

##      ICC1 ICC1_lower ICC1_upper  ICC3  ICC3_lower ICC3_upper ICC3k ICC3k_lower ICC3k_upper               var
## 1  -0.011 -0.4163157  0.4013330 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192         -0.5 drug
## 2  -0.233 -0.5842313  0.1969920 0.048 -0.37265010  0.4522652 0.092 -1.18801360   0.6228410         -1.5 drug
## 3   0.180 -0.2452362  0.5498219 0.180 -0.25226520  0.5519705 0.305 -0.67474510   0.7113157          0.5 drug
## 4   0.000 -0.4070147  0.4106858 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192          1.5 drug
## 5  -0.042 -0.4417737  0.3748305 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192              drug
## 6   0.281 -0.1420422  0.6202767 0.281 -0.14937578  0.6221709 0.439 -0.35121451   0.7670843      -0.5 neutral
## 7   0.128 -0.2940845  0.5119968 0.129 -0.30031985  0.5147496 0.228 -0.85844898   0.6796497      -1.5 neutral
## 8   0.000 -0.4070147  0.4106858 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192       0.5 neutral
## 9   0.000 -0.4070147  0.4106858 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192       1.5 neutral
## 10  0.000 -0.4070147  0.4106858 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192           neutral
## 11  0.000 -0.4070147  0.4106858 0.000 -0.41324703  0.4132470 0.000 -1.40858948   0.5848192  -0.5DrugvNeutral
## 12  0.050 -0.3645727  0.4512863 0.440  0.03306215  0.7221461 0.611  0.06400805   0.8386583  -1.5DrugvNeutral
## 13  0.165 -0.2592465  0.5392968 0.165 -0.26622191  0.5414807 0.284 -0.72561967   0.7025462   0.5DrugvNeutral
## 14  0.202 -0.2236710  0.5655438 0.202 -0.23077683  0.5676383 0.336 -0.60002570   0.7241955   1.5DrugvNeutral
## 15 -0.131 -0.5104297  0.2960213 0.125 -0.30351702  0.5121592 0.223 -0.87157053   0.6773879 DrugvNeutralEarly
## 16  0.236 -0.1889639  0.5896970 0.236 -0.19617868  0.5917054 0.382 -0.48811515   0.7434861  DrugvNeutralLate
## 17  0.218 -0.2078345  0.5767352 0.317 -0.11044641  0.6458480 0.482 -0.24831874   0.7848210      DrugvNeutral
## Saving 7 x 5 in image

## [1] "###tDCS MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 168.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.4842 -0.5292  0.0120  0.5706  4.7914 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.00000  0.0000  
##  id           (Intercept) 0.01270  0.1127  
##  Residual                 0.06946  0.2635  
## Number of obs: 520, groups:  condition:id, 130; id, 65
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)         -0.02950    0.02151 123.63325  -1.372 0.172678    
## conditiondrug        0.11960    0.02311 452.00000   5.174 3.45e-07 ***
## time                 0.02737    0.01462 452.00000   1.872 0.061827 .  
## conditiondrug:time  -0.07181    0.02067 452.00000  -3.473 0.000564 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.537              
## time         0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707
## [1] "###Neurocaps MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 140.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.6720 -0.4571 -0.0145  0.4209  6.4603 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.00000  0.0000  
##  id           (Intercept) 0.02804  0.1674  
##  Residual                 0.08492  0.2914  
## Number of obs: 232, groups:  condition:id, 58; id, 29
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)  
## (Intercept)         -0.04553    0.04122  45.01645  -1.105   0.2752  
## conditiondrug        0.08482    0.03826 200.00000   2.217   0.0278 *
## time                 0.05273    0.02420 200.00000   2.179   0.0305 *
## conditiondrug:time  -0.07873    0.03422 200.00000  -2.300   0.0225 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.464              
## time         0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707
## [1] "###Neurocaps OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 156.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.5623 -0.5930  0.0633  0.5813  2.9577 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.00000  0.0000  
##  id           (Intercept) 0.02535  0.1592  
##  Residual                 0.11549  0.3398  
## Number of obs: 176, groups:  condition:id, 44; id, 22
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)   
## (Intercept)         -0.01949    0.04965  38.30370  -0.393  0.69678   
## conditiondrug        0.10844    0.05123 151.00000   2.117  0.03592 * 
## time                 0.03162    0.03240 151.00000   0.976  0.33072   
## conditiondrug:time  -0.12112    0.04582 151.00000  -2.643  0.00908 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.516              
## time         0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707
## [1] "###tACS OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 75.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.3871 -0.5269  0.0218  0.5408  2.3326 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.00000  0.0000  
##  id           (Intercept) 0.01807  0.1344  
##  Residual                 0.07131  0.2670  
## Number of obs: 176, groups:  condition:id, 44; id, 22
## 
## Fixed effects:
##                      Estimate Std. Error         df t value Pr(>|t|)
## (Intercept)         -0.015030   0.040394  36.609319  -0.372    0.712
## conditiondrug       -0.006365   0.040257 151.000000  -0.158    0.875
## time                -0.011709   0.025461 151.000000  -0.460    0.646
## conditiondrug:time   0.027189   0.036007 151.000000   0.755    0.451
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.498              
## time         0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707

RVStriatum

knitr::include_graphics("RVStri_zm1.png")

plot_combination('6', 'RVSTRI', limits_y = c(-0.2, 0.3))
## [1] "###Condition z-scores###"
##                    tdcs neurocaps_mcr neurocaps_ocr      tacs
## tdcs                 NA     0.5332700     0.7414731 1.9582658
## neurocaps_mcr 0.5332700            NA     0.2394436 1.2620420
## neurocaps_ocr 0.7414731     0.2394436            NA 0.9295674
## tacs          1.9582658     1.2620420     0.9295674        NA
## [1] "###Condition by time z-scores###"
##                     tdcs neurocaps_mcr neurocaps_ocr      tacs
## tdcs                  NA    -0.7847964     0.7936411 -2.304104
## neurocaps_mcr -0.7847964            NA     1.2817258 -1.361958
## neurocaps_ocr  0.7936411     1.2817258            NA -2.460801
## tacs          -2.3041039    -1.3619583    -2.4608006        NA
## [1] "###Condition p-values###"
##                     tdcs neurocaps_mcr neurocaps_ocr       tacs
## tdcs                  NA     0.5938467     0.4584066 0.05019883
## neurocaps_mcr 0.59384666            NA     0.8107616 0.20693369
## neurocaps_ocr 0.45840664     0.8107616            NA 0.35259512
## tacs          0.05019883     0.2069337     0.3525951         NA
## [1] "###Condition by time p-values###"
##                     tdcs neurocaps_mcr neurocaps_ocr       tacs
## tdcs                  NA     0.4325729    0.42740442 0.02121681
## neurocaps_mcr 0.43257295            NA    0.19993887 0.17321106
## neurocaps_ocr 0.42740442     0.1999389            NA 0.01386274
## tacs          0.02121681     0.1732111    0.01386274         NA
## [1] "one"
## [1] "one"
## [1] "one"
## [1] "one"
## Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

##      ICC1  ICC1_lower ICC1_upper  ICC3  ICC3_lower ICC3_upper ICC3k ICC3k_lower ICC3k_upper               var
## 1  -0.013 -0.41814081  0.3994775 0.000 -0.41324703  0.4132470 0.000  -1.4085895   0.5848192         -0.5 drug
## 2  -0.183 -0.54883570  0.2465632 0.000 -0.41324703  0.4132470 0.000  -1.4085895   0.5848192         -1.5 drug
## 3   0.173 -0.25189290  0.5448524 0.173 -0.25889682  0.5470178 0.295  -0.6986795   0.7071900          0.5 drug
## 4   0.209 -0.21594380  0.5710411 0.209 -0.22307541  0.5731163 0.346  -0.5742524   0.7286381          1.5 drug
## 5  -0.004 -0.41043615  0.4072652 0.000 -0.41324703  0.4132470 0.000  -1.4085895   0.5848192              drug
## 6   0.466  0.07263964  0.7357876 0.466  0.06518258  0.7371991 0.636   0.1223876   0.8487215      -0.5 neutral
## 7   0.183 -0.24181183  0.5523566 0.183 -0.24885350  0.5544966 0.310  -0.6625964   0.7134098      -1.5 neutral
## 8   0.000 -0.40701472  0.4106858 0.000 -0.41324703  0.4132470 0.000  -1.4085895   0.5848192       0.5 neutral
## 9   0.000 -0.40701472  0.4106858 0.000 -0.41324703  0.4132470 0.000  -1.4085895   0.5848192       1.5 neutral
## 10  0.000 -0.40701472  0.4106858 0.000 -0.41324703  0.4132470 0.000  -1.4085895   0.5848192           neutral
## 11 -0.007 -0.41276396  0.4049248 0.000 -0.41324703  0.4132470 0.000  -1.4085895   0.5848192  -0.5DrugvNeutral
## 12 -0.103 -0.48941402  0.3213702 0.030 -0.38776739  0.4380956 0.059  -1.2667322   0.6092719  -1.5DrugvNeutral
## 13  0.204 -0.22138222  0.5671794 0.204 -0.22849574  0.5692682 0.339  -0.5923383   0.7255206   0.5DrugvNeutral
## 14  0.081 -0.33694092  0.4760205 0.093 -0.33326886  0.4873205 0.170  -0.9997099   0.6552999   1.5DrugvNeutral
## 15 -0.089 -0.47856170  0.3340204 0.046 -0.37437359  0.4506705 0.088  -1.1967960   0.6213272 DrugvNeutralEarly
## 16  0.149 -0.27471736  0.5273771 0.149 -0.28163020  0.5296003 0.259  -0.7840814   0.6924689  DrugvNeutralLate
## 17  0.308 -0.11316548  0.6380107 0.319 -0.10837702  0.6470672 0.484  -0.2431006   0.7857205      DrugvNeutral
## Saving 7 x 5 in image

## [1] "###tDCS MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 91.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.1701 -0.4820  0.0024  0.5321  4.7283 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.00000  0.0000  
##  id           (Intercept) 0.01108  0.1053  
##  Residual                 0.05978  0.2445  
## Number of obs: 520, groups:  condition:id, 130; id, 65
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)         -0.04614    0.02001 123.11628  -2.306 0.022788 *  
## conditiondrug        0.10634    0.02144 452.00000   4.959    1e-06 ***
## time                 0.01734    0.01356 452.00000   1.279 0.201721    
## conditiondrug:time  -0.06362    0.01918 452.00000  -3.317 0.000984 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.536              
## time         0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707
## [1] "###Neurocaps MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 77.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6241 -0.5212 -0.0686  0.5373  5.1892 
## 
## Random effects:
##  Groups       Name        Variance  Std.Dev. 
##  condition:id (Intercept) 6.361e-18 2.522e-09
##  id           (Intercept) 1.236e-02 1.112e-01
##  Residual                 6.758e-02 2.600e-01
## Number of obs: 232, groups:  condition:id, 58; id, 29
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)  
## (Intercept)         -0.07995    0.03176  54.10745  -2.517   0.0148 *
## conditiondrug        0.08485    0.03413 200.00000   2.486   0.0137 *
## time                 0.01649    0.02159 200.00000   0.764   0.4458  
## conditiondrug:time  -0.03532    0.03053 200.00000  -1.157   0.2487  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.537              
## time         0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707
## [1] "###Neurocaps OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 78.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3397 -0.5610  0.1025  0.5943  2.4400 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.00000  0.0000  
##  id           (Intercept) 0.01485  0.1219  
##  Residual                 0.07369  0.2715  
## Number of obs: 176, groups:  condition:id, 44; id, 22
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)   
## (Intercept)         -0.02814    0.03889  39.35223  -0.724  0.47361   
## conditiondrug        0.07209    0.04092 151.00000   1.762  0.08017 . 
## time                 0.00827    0.02588 151.00000   0.320  0.74977   
## conditiondrug:time  -0.09641    0.03660 151.00000  -2.634  0.00932 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.526              
## time         0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707
## [1] "###tACS OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 64.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.2814 -0.4684  0.0312  0.5530  2.7120 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.00000  0.0000  
##  id           (Intercept) 0.02035  0.1426  
##  Residual                 0.06581  0.2565  
## Number of obs: 176, groups:  condition:id, 44; id, 22
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)
## (Intercept)         -0.04331    0.04090  34.43523  -1.059    0.297
## conditiondrug        0.01975    0.03867 151.00000   0.511    0.610
## time                -0.02562    0.02446 151.00000  -1.047    0.297
## conditiondrug:time   0.02752    0.03459 151.00000   0.796    0.428
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.473              
## time         0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707

RAmy

knitr::include_graphics("RAmy_zm16.png")

plot_combination('8', 'RAmy', limits_y = c(-0.05, 0.75))
## [1] "###Condition z-scores###"
##                    tdcs neurocaps_mcr neurocaps_ocr      tacs
## tdcs                 NA     0.7234911     2.4657736 2.7482396
## neurocaps_mcr 0.7234911            NA     1.3756300 1.8963280
## neurocaps_ocr 2.4657736     1.3756300            NA 0.7971716
## tacs          2.7482396     1.8963280     0.7971716        NA
## [1] "###Condition by time z-scores###"
##                     tdcs neurocaps_mcr neurocaps_ocr       tacs
## tdcs                  NA    -0.6276469     0.3958587 -1.4661307
## neurocaps_mcr -0.6276469            NA     0.8888891 -0.8061678
## neurocaps_ocr  0.3958587     0.8888891            NA -1.6096743
## tacs          -1.4661307    -0.8061678    -1.6096743         NA
## [1] "###Condition p-values###"
##                      tdcs neurocaps_mcr neurocaps_ocr        tacs
## tdcs                   NA    0.46937823    0.01367177 0.005991621
## neurocaps_mcr 0.469378227            NA    0.16893623 0.057916683
## neurocaps_ocr 0.013671771    0.16893623            NA 0.425351380
## tacs          0.005991621    0.05791668    0.42535138          NA
## [1] "###Condition by time p-values###"
##                    tdcs neurocaps_mcr neurocaps_ocr      tacs
## tdcs                 NA     0.5302353     0.6922093 0.1426127
## neurocaps_mcr 0.5302353            NA     0.3740627 0.4201461
## neurocaps_ocr 0.6922093     0.3740627            NA 0.1074690
## tacs          0.1426127     0.4201461     0.1074690        NA
## [1] "one"
## [1] "one"
## [1] "one"
## [1] "one"
## Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

##      ICC1 ICC1_lower ICC1_upper  ICC3 ICC3_lower ICC3_upper ICC3k ICC3k_lower ICC3k_upper               var
## 1   0.000 -0.4070147  0.4106858 0.000 -0.4132470  0.4132470 0.000  -1.4085895   0.5848192         -0.5 drug
## 2  -0.028 -0.4299547  0.3873026 0.001 -0.4123668  0.4141265 0.002  -1.4034836   0.5856994         -1.5 drug
## 3   0.080 -0.3379906  0.4751038 0.081 -0.3436572  0.4783232 0.150  -1.0471881   0.6471158          0.5 drug
## 4   0.029 -0.3828382  0.4342177 0.029 -0.3892143  0.4367175 0.056  -1.2744711   0.6079379          1.5 drug
## 5   0.000 -0.4070147  0.4106858 0.000 -0.4132470  0.4132470 0.000  -1.4085895   0.5848192              drug
## 6  -0.169 -0.5392473  0.2593116 0.000 -0.4132470  0.4132470 0.000  -1.4085895   0.5848192      -0.5 neutral
## 7   0.213 -0.2123710  0.5735592 0.213 -0.2195142  0.5756255 0.351  -0.5625066   0.7306628      -1.5 neutral
## 8   0.221 -0.2048928  0.5787821 0.227 -0.2057729  0.5851693 0.370  -0.5181713   0.7383051       0.5 neutral
## 9   0.000 -0.4070147  0.4106858 0.000 -0.4132470  0.4132470 0.000  -1.4085895   0.5848192       1.5 neutral
## 10 -0.023 -0.4256577  0.3917642 0.000 -0.4132470  0.4132470 0.000  -1.4085895   0.5848192           neutral
## 11  0.000 -0.4070147  0.4106858 0.000 -0.4132470  0.4132470 0.000  -1.4085895   0.5848192  -0.5DrugvNeutral
## 12  0.169 -0.2555671  0.5420853 0.214 -0.2181521  0.5765813 0.353  -0.5580423   0.7314324  -1.5DrugvNeutral
## 13  0.295 -0.1275845  0.6292541 0.295 -0.1349480  0.6311138 0.455  -0.3119998   0.7738440   0.5DrugvNeutral
## 14  0.000 -0.4070147  0.4106858 0.000 -0.4132470  0.4132470 0.000  -1.4085895   0.5848192   1.5DrugvNeutral
## 15  0.244 -0.1812008  0.5949153 0.324 -0.1030006  0.6502167 0.489  -0.2296560   0.7880380 DrugvNeutralEarly
## 16  0.137 -0.2863589  0.5181947 0.137 -0.2932223  0.5204475 0.240  -0.8297439   0.6845978  DrugvNeutralLate
## 17  0.144 -0.2793620  0.5237359 0.144 -0.2862554  0.5259709 0.252  -0.8021227   0.6893590      DrugvNeutral
## Saving 7 x 5 in image

## [1] "###tDCS MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 598.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.4349 -0.5583  0.0212  0.5374  5.3478 
## 
## Random effects:
##  Groups       Name        Variance  Std.Dev.
##  condition:id (Intercept) 0.0007294 0.02701 
##  id           (Intercept) 0.0467183 0.21614 
##  Residual                 0.1527508 0.39083 
## Number of obs: 520, groups:  condition:id, 130; id, 65
## 
## Fixed effects:
##                      Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)          0.165024   0.036297  98.641800   4.547 1.55e-05 ***
## conditiondrug        0.292493   0.034604  63.999973   8.453 5.15e-12 ***
## time                -0.004978   0.021680 387.999967  -0.230 0.818492    
## conditiondrug:time  -0.105186   0.030659 387.999967  -3.431 0.000667 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.477              
## time         0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707
## [1] "###Neurocaps MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 220.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6456 -0.6269 -0.0239  0.5272  3.8607 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.01345  0.1160  
##  id           (Intercept) 0.01775  0.1332  
##  Residual                 0.11989  0.3462  
## Number of obs: 232, groups:  condition:id, 58; id, 29
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)          0.16800    0.04593  51.65200   3.658 0.000597 ***
## conditiondrug        0.24565    0.05472  28.00000   4.489 0.000112 ***
## time                 0.04765    0.02875 172.00000   1.657 0.099304 .  
## conditiondrug:time  -0.07322    0.04066 172.00000  -1.801 0.073518 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.596              
## time         0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707
## [1] "###Neurocaps OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 144.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2559 -0.5700  0.0244  0.5475  3.6217 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.0000   0.0000  
##  id           (Intercept) 0.0294   0.1715  
##  Residual                 0.1056   0.3249  
## Number of obs: 176, groups:  condition:id, 44; id, 22
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)   
## (Intercept)          0.15363    0.05036  35.55075   3.051  0.00430 **
## conditiondrug        0.14462    0.04898 151.00000   2.953  0.00366 **
## time                 0.01792    0.03098 151.00000   0.578  0.56386   
## conditiondrug:time  -0.12635    0.04381 151.00000  -2.884  0.00450 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.486              
## time         0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707
## [1] "###tACS OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time + (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 179.1
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.90348 -0.50601 -0.00016  0.53829  2.35498 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.02513  0.1585  
##  id           (Intercept) 0.02141  0.1463  
##  Residual                 0.11978  0.3461  
## Number of obs: 176, groups:  condition:id, 44; id, 22
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)          0.24023    0.05897  38.94810   4.074  0.00022 ***
## conditiondrug        0.07601    0.07076  21.00000   1.074  0.29491    
## time                -0.04487    0.03300 130.00000  -1.360  0.17622    
## conditiondrug:time  -0.02332    0.04667 130.00000  -0.500  0.61813    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.600              
## time         0.000  0.000       
## cndtndrg:tm  0.000  0.000 -0.707

Self-Report and RTs

Craving

plot_combination('craving', 'Craving', limits_y = c(1, 4), imaging = FALSE, forest_range = c(-0.5, 2.5))
## [1] "###Condition z-scores###"
##                     tdcs neurocaps_mcr neurocaps_ocr       tacs
## tdcs                  NA     0.1576549    -3.5832372 -3.6285460
## neurocaps_mcr  0.1576549            NA    -2.8131019 -2.8858501
## neurocaps_ocr -3.5832372    -2.8131019            NA -0.1371176
## tacs          -3.6285460    -2.8858501    -0.1371176         NA
## [1] "###Condition by time z-scores###"
##                     tdcs neurocaps_mcr neurocaps_ocr        tacs
## tdcs                  NA   -1.36784194    -0.6971554 -1.10796210
## neurocaps_mcr -1.3678419            NA     0.4683281  0.08924022
## neurocaps_ocr -0.6971554    0.46832809            NA -0.34567535
## tacs          -1.1079621    0.08924022    -0.3456754          NA
## [1] "###Condition p-values###"
##                       tdcs neurocaps_mcr neurocaps_ocr        tacs
## tdcs                    NA    0.87472870  0.0003393621 0.000285022
## neurocaps_mcr 0.8747287022            NA  0.0049066102 0.003903580
## neurocaps_ocr 0.0003393621    0.00490661            NA 0.890937842
## tacs          0.0002850220    0.00390358  0.8909378416          NA
## [1] "###Condition by time p-values###"
##                    tdcs neurocaps_mcr neurocaps_ocr      tacs
## tdcs                 NA     0.1713616     0.4857055 0.2678782
## neurocaps_mcr 0.1713616            NA     0.6395500 0.9288910
## neurocaps_ocr 0.4857055     0.6395500            NA 0.7295867
## tacs          0.2678782     0.9288910     0.7295867        NA
## [1] "one"
## [1] "one"
## [1] "one"
## [1] "one"
## Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

##     ICC1  ICC1_lower ICC1_upper  ICC3   ICC3_lower ICC3_upper ICC3k ICC3k_lower ICC3k_upper               var
## 1  0.640  0.31473935  0.8318957 0.648  0.321032941  0.8372305 0.787  0.48603321   0.9114049         -0.5 drug
## 2  0.210 -0.21495602  0.5717387 0.210 -0.222090866  0.5738115 0.348 -0.57099437   0.7291998         -1.5 drug
## 3  0.290 -0.13220264  0.6264082 0.290 -0.139556940  0.6282790 0.450 -0.32438391   0.7717092          0.5 drug
## 4  0.408  0.00172231  0.7014687 0.408 -0.005770302  0.7030322 0.580 -0.01160758   0.8256241          1.5 drug
## 5  0.524  0.14929176  0.7694099 0.524  0.141958008  0.7706655 0.688  0.24862212   0.8704812              drug
## 6  0.593  0.24556160  0.8072723 0.601  0.250377263  0.8126749 0.751  0.40048275   0.8966582      -0.5 neutral
## 7  0.457  0.06109960  0.7304258 0.478  0.081211991  0.7444697 0.647  0.15022399   0.8535198      -1.5 neutral
## 8  0.256 -0.16844405  0.6033506 0.256 -0.175714827  0.6053087 0.408 -0.42634475   0.7541337       0.5 neutral
## 9  0.227 -0.19834601  0.5833024 0.227 -0.205533104  0.5853339 0.370 -0.51741137   0.7384361       1.5 neutral
## 10 0.471  0.07915571  0.7387789 0.471  0.071705701  0.7401768 0.640  0.13381603   0.8506915           neutral
## 11 0.549  0.18239307  0.7829393 0.599  0.246868541  0.8114011 0.749  0.39598167   0.8958823  -0.5DrugvNeutral
## 12 0.169 -0.25556597  0.5420861 0.169 -0.262555756  0.5442607 0.289 -0.71206944   0.7048819  -1.5DrugvNeutral
## 13 0.450  0.05319355  0.7267043 0.456  0.052062876  0.7311302 0.626  0.09897294   0.8446854   0.5DrugvNeutral
## 14 0.585  0.23387444  0.8029099 0.585  0.226779286  0.8040033 0.738  0.36971489   0.8913546   1.5DrugvNeutral
## 15 0.470  0.07798528  0.7382435 0.470  0.070533954  0.7396438 0.640  0.13177341   0.8503394 DrugvNeutralEarly
## 16 0.747  0.48816843  0.8858083 0.753  0.493637658  0.8895781 0.859  0.66098716   0.9415626  DrugvNeutralLate
## 17 0.698  0.40666757  0.8617427 0.698  0.400395025  0.8625346 0.822  0.57183154   0.9261944      DrugvNeutral
## Saving 7 x 5 in image

## [1] "###tDCS MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time +  (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 1014
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.06640 -0.53291  0.05142  0.55519  3.06609 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.1242   0.3524  
##  id           (Intercept) 0.3488   0.5906  
##  Residual                 0.2824   0.5314  
## Number of obs: 494, groups:  condition:id, 129; id, 65
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)          1.69379    0.09207  91.37980   18.40  < 2e-16 ***
## conditiondrug        1.31929    0.07862  63.39930   16.78  < 2e-16 ***
## time                 0.18236    0.03106 368.81187    5.87  9.7e-09 ***
## conditiondrug:time  -0.10600    0.04327 367.10038   -2.45   0.0148 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.427              
## time        -0.013  0.013       
## cndtndrg:tm  0.009 -0.015 -0.718
## [1] "###Neurocaps MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time +  (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 370.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3459 -0.4088 -0.0669  0.1513  3.7800 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.36796  0.6066  
##  id           (Intercept) 0.07349  0.2711  
##  Residual                 0.14943  0.3866  
## Number of obs: 228, groups:  condition:id, 58; id, 29
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)          1.08915    0.12866  54.96792   8.465 1.54e-11 ***
## conditiondrug        1.29014    0.16737  27.95315   7.708 2.16e-08 ***
## time                 0.04062    0.03221 168.23056   1.261    0.209    
## conditiondrug:time  -0.01991    0.04571 168.28819  -0.436    0.664    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.651              
## time        -0.001  0.001       
## cndtndrg:tm  0.001 -0.003 -0.705
## [1] "###Neurocaps OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time +  (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 296.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.8107 -0.5292  0.0609  0.4771  2.8296 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.1905   0.4365  
##  id           (Intercept) 0.1031   0.3212  
##  Residual                 0.1874   0.4329  
## Number of obs: 174, groups:  condition:id, 44; id, 22
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)          1.38243    0.12454  38.34890  11.101 1.53e-13 ***
## conditiondrug        1.91698    0.14711  21.01962  13.031 1.55e-11 ***
## time                 0.14109    0.04182 128.19406   3.374 0.000981 ***
## conditiondrug:time  -0.05492    0.05914 128.18996  -0.929 0.354867    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.591              
## time        -0.007  0.006       
## cndtndrg:tm  0.005 -0.009 -0.707
## [1] "###tACS OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time +  (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 291.3
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.64887 -0.62209  0.04536  0.61595  2.45138 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.2149   0.4636  
##  id           (Intercept) 0.1002   0.3165  
##  Residual                 0.1777   0.4215  
## Number of obs: 173, groups:  condition:id, 44; id, 22
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)          1.38609    0.12795  39.04699  10.833 2.49e-13 ***
## conditiondrug        1.94616    0.15383  20.95862  12.651 2.81e-11 ***
## time                 0.17284    0.04025 127.16548   4.295 3.44e-05 ***
## conditiondrug:time  -0.02645    0.05730 127.31114  -0.462    0.645    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.600              
## time        -0.002  0.002       
## cndtndrg:tm  0.002 -0.004 -0.702

Craving RT

plot_combination('craving_rt', 'CravingRT', limits_y = c(1, 3.5), imaging = FALSE)
## [1] "###Condition z-scores###"
##                     tdcs neurocaps_mcr neurocaps_ocr       tacs
## tdcs                  NA     -2.385365      1.594898 -0.8158398
## neurocaps_mcr -2.3853647            NA      3.496007  1.4708419
## neurocaps_ocr  1.5948985      3.496007            NA -2.1848318
## tacs          -0.8158398      1.470842     -2.184832         NA
## [1] "###Condition by time z-scores###"
##                      tdcs neurocaps_mcr neurocaps_ocr      tacs
## tdcs                   NA    0.09861327    -0.1295825 0.3083629
## neurocaps_mcr  0.09861327            NA    -0.1923227 0.1838783
## neurocaps_ocr -0.12958247   -0.19232270            NA 0.3477431
## tacs           0.30836289    0.18387829     0.3477431        NA
## [1] "###Condition p-values###"
##                    tdcs neurocaps_mcr neurocaps_ocr       tacs
## tdcs                 NA  0.0170621969  0.1107349440 0.41459174
## neurocaps_mcr 0.0170622            NA  0.0004722765 0.14133389
## neurocaps_ocr 0.1107349  0.0004722765            NA 0.02890118
## tacs          0.4145917  0.1413338910  0.0289011789         NA
## [1] "###Condition by time p-values###"
##                    tdcs neurocaps_mcr neurocaps_ocr      tacs
## tdcs                 NA     0.9214453     0.8968968 0.7578062
## neurocaps_mcr 0.9214453            NA     0.8474894 0.8541089
## neurocaps_ocr 0.8968968     0.8474894            NA 0.7280331
## tacs          0.7578062     0.8541089     0.7280331        NA
## [1] "one"
## [1] "one"
## [1] "one"
## [1] "one"
## Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

##      ICC1  ICC1_lower ICC1_upper  ICC3   ICC3_lower ICC3_upper ICC3k ICC3k_lower ICC3k_upper               var
## 1   0.010 -0.39852452  0.4190756 0.010 -0.404808316  0.4216154 0.020 -1.36026200   0.5931497         -0.5 drug
## 2  -0.037 -0.43731361  0.3795720 0.000 -0.413247031  0.4132470 0.000 -1.40858948   0.5848192         -1.5 drug
## 3   0.000 -0.40701472  0.4106858 0.000 -0.413247031  0.4132470 0.000 -1.40858948   0.5848192          0.5 drug
## 4   0.439  0.03843318  0.7196487 0.439  0.030949557  0.7211326 0.610  0.06004088   0.8379745          1.5 drug
## 5   0.036 -0.37652310  0.4401861 0.036 -0.382935331  0.4426698 0.070 -1.24115137   0.6136814              drug
## 6  -0.012 -0.41685610  0.4007843 0.036 -0.383271168  0.4423533 0.069 -1.24291633   0.6133772      -0.5 neutral
## 7  -0.140 -0.51762725  0.2870705 0.064 -0.358693986  0.4649868 0.120 -1.11863596   0.6348000      -1.5 neutral
## 8   0.471  0.07975234  0.7390515 0.471  0.072303002  0.7404481 0.641  0.13485554   0.8508707       0.5 neutral
## 9   0.325 -0.09494472  0.6488050 0.359 -0.063215014  0.6727456 0.529 -0.13496163   0.8043609       1.5 neutral
## 10  0.355 -0.06111637  0.6680822 0.564  0.196853469  0.7926416 0.721  0.32895166   0.8843280           neutral
## 11  0.021 -0.38937000  0.4279682 0.033 -0.385178854  0.4405514 0.065 -1.25297855   0.6116427  -0.5DrugvNeutral
## 12 -0.012 -0.41669523  0.4009477 0.000 -0.413247031  0.4132470 0.000 -1.40858948   0.5848192  -1.5DrugvNeutral
## 13  0.188 -0.23692188  0.5559509 0.188 -0.243981312  0.5580785 0.317 -0.64543725   0.7163677   0.5DrugvNeutral
## 14  0.372 -0.04139961  0.6788837 0.420  0.007584464  0.7097232 0.591  0.01505475   0.8302200   1.5DrugvNeutral
## 15 -0.055 -0.45212917  0.3636542 0.000 -0.413247031  0.4132470 0.000 -1.40858948   0.5848192 DrugvNeutralEarly
## 16  0.116 -0.30521874  0.5029143 0.125 -0.303592532  0.5120978 0.223 -0.87188189   0.6773342  DrugvNeutralLate
## 17 -0.054 -0.45070979  0.3652001 0.000 -0.413247031  0.4132470 0.000 -1.40858948   0.5848192      DrugvNeutral
## Saving 7 x 5 in image

## [1] "###tDCS MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time +  (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 1320.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1302 -0.6149 -0.2205  0.4779  3.4277 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.1413   0.3759  
##  id           (Intercept) 0.1186   0.3444  
##  Residual                 0.6601   0.8125  
## Number of obs: 494, groups:  condition:id, 129; id, 65
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)          2.09015    0.08225 119.66236  25.411  < 2e-16 ***
## conditiondrug       -0.24541    0.09890  64.12233  -2.481   0.0157 *  
## time                -0.22534    0.04734 373.25097  -4.760 2.78e-06 ***
## conditiondrug:time   0.09599    0.06602 370.24830   1.454   0.1468    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.606              
## time        -0.019  0.015       
## cndtndrg:tm  0.014 -0.018 -0.717
## [1] "###Neurocaps MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time +  (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 570.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8201 -0.5548 -0.2002  0.3057  4.1821 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.06786  0.2605  
##  id           (Intercept) 0.14683  0.3832  
##  Residual                 0.55775  0.7468  
## Number of obs: 228, groups:  condition:id, 58; id, 29
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)          1.52463    0.11116  48.04003  13.716  < 2e-16 ***
## conditiondrug        0.12622    0.12038  27.94125   1.049    0.303    
## time                -0.29659    0.06221 168.51298  -4.768 4.01e-06 ***
## conditiondrug:time   0.08512    0.08824 168.78397   0.965    0.336    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.545              
## time        -0.003  0.003       
## cndtndrg:tm  0.002 -0.007 -0.705
## [1] "###Neurocaps OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time +  (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 481.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4627 -0.6814 -0.1119  0.5459  2.8380 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.02094  0.1447  
##  id           (Intercept) 0.19729  0.4442  
##  Residual                 0.76429  0.8742  
## Number of obs: 174, groups:  condition:id, 44; id, 22
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)          2.23722    0.13681  34.31196  16.352  < 2e-16 ***
## conditiondrug       -0.51833    0.13965  20.74077  -3.712  0.00131 ** 
## time                -0.44598    0.08432 128.93920  -5.289 5.12e-07 ***
## conditiondrug:time   0.11365    0.11926 128.92380   0.953  0.34237    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.510              
## time        -0.012  0.012       
## cndtndrg:tm  0.008 -0.016 -0.707
## [1] "###tACS OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time +  (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 391.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6610 -0.6255 -0.2142  0.4807  2.8207 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.04606  0.2146  
##  id           (Intercept) 0.23621  0.4860  
##  Residual                 0.40815  0.6389  
## Number of obs: 173, groups:  condition:id, 44; id, 22
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)          1.75887    0.13240  30.35512  13.284 3.47e-14 ***
## conditiondrug       -0.12052    0.11684  21.27910  -1.032  0.31389    
## time                -0.20351    0.06099 127.29143  -3.337  0.00111 ** 
## conditiondrug:time   0.06236    0.08679 127.67725   0.719  0.47372    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.439              
## time        -0.003  0.004       
## cndtndrg:tm  0.002 -0.008 -0.703

Box RT

plot_combination('box_rt', 'CravingRT', limits_y = c(0.5, 3.5), imaging = FALSE)
## [1] "###Condition z-scores###"
##                    tdcs neurocaps_mcr neurocaps_ocr       tacs
## tdcs                 NA     -0.789830     2.0756962  1.6805754
## neurocaps_mcr -0.789830            NA     2.4623103  2.1387617
## neurocaps_ocr  2.075696      2.462310            NA -0.6546852
## tacs           1.680575      2.138762    -0.6546852         NA
## [1] "###Condition by time z-scores###"
##                     tdcs neurocaps_mcr neurocaps_ocr       tacs
## tdcs                  NA      2.534853      1.091190 -0.1697659
## neurocaps_mcr  2.5348532            NA     -1.416087 -2.4481921
## neurocaps_ocr  1.0911901     -1.416087            NA -1.1194922
## tacs          -0.1697659     -2.448192     -1.119492         NA
## [1] "###Condition p-values###"
##                     tdcs neurocaps_mcr neurocaps_ocr       tacs
## tdcs                  NA    0.42962703    0.03792206 0.09284541
## neurocaps_mcr 0.42962703            NA    0.01380452 0.03245497
## neurocaps_ocr 0.03792206    0.01380452            NA 0.51267043
## tacs          0.09284541    0.03245497    0.51267043         NA
## [1] "###Condition by time p-values###"
##                     tdcs neurocaps_mcr neurocaps_ocr       tacs
## tdcs                  NA    0.01124945     0.2751893 0.86519428
## neurocaps_mcr 0.01124945            NA     0.1567501 0.01435751
## neurocaps_ocr 0.27518926    0.15675014            NA 0.26293022
## tacs          0.86519428    0.01435751     0.2629302         NA
## [1] "one"
## [1] "one"
## [1] "one"
## [1] "one"
## Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

##      ICC1  ICC1_lower ICC1_upper  ICC3  ICC3_lower ICC3_upper ICC3k ICC3k_lower ICC3k_upper               var
## 1   0.000 -0.40701472  0.4106858 0.000 -0.41324703  0.4132470 0.000  -1.4085895   0.5848192         -0.5 drug
## 2   0.671  0.36328604  0.8480396 0.671  0.35676460  0.8489036 0.803   0.5259049   0.9182778         -1.5 drug
## 3   0.195 -0.23004587  0.5609553 0.292 -0.13806489  0.6291989 0.452  -0.3203603   0.7724028          0.5 drug
## 4   0.045 -0.36852838  0.4476398 0.049 -0.37193485  0.4529254 0.093  -1.1843830   0.6234668          1.5 drug
## 5   0.498  0.11444423  0.7545390 0.498  0.10704348  0.7558644 0.665   0.1933862   0.8609599              drug
## 6   0.276 -0.14787226  0.6165991 0.276 -0.15446338  0.6189687 0.433  -0.3653618   0.7646457      -0.5 neutral
## 7   0.075 -0.34282183  0.4708619 0.075 -0.34941686  0.4732592 0.139  -1.0741651   0.6424657      -1.5 neutral
## 8   0.810  0.60106781  0.9159504 0.810  0.59626054  0.9164458 0.895   0.7470717   0.9564015       0.5 neutral
## 9   0.249 -0.17599277  0.5983798 0.306 -0.12247901  0.6386818 0.469  -0.2791478   0.7795068       1.5 neutral
## 10  0.365 -0.04945140  0.6745100 0.365 -0.05692284  0.6761878 0.535  -0.1207173   0.8068163           neutral
## 11  0.000 -0.40701472  0.4106858 0.000 -0.41324703  0.4132470 0.000  -1.4085895   0.5848192  -0.5DrugvNeutral
## 12  0.715  0.43418603  0.8701040 0.715  0.42808612  0.8708513 0.834   0.5995242   0.9309680  -1.5DrugvNeutral
## 13 -0.114 -0.49824159  0.3108622 0.000 -0.41324703  0.4132470 0.000  -1.4085895   0.5848192   0.5DrugvNeutral
## 14  0.325 -0.09457070  0.6490235 0.457  0.05410178  0.7320805 0.627   0.1026500   0.8453192   1.5DrugvNeutral
## 15  0.224 -0.20099234  0.5814810 0.227 -0.20541124  0.5854176 0.370  -0.5170253   0.7385027 DrugvNeutralEarly
## 16  0.467  0.07401810  0.7364226 0.467  0.06656248  0.7378312 0.637   0.1248168   0.8491403  DrugvNeutralLate
## 17  0.349 -0.06786693  0.6643120 0.349 -0.07532117  0.6660319 0.517  -0.1629131   0.7995428      DrugvNeutral
## Saving 7 x 5 in image

## [1] "###tDCS MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time +  (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 1086.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8130 -0.5455 -0.1869  0.2686  3.4076 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.1706   0.4130  
##  id           (Intercept) 0.1763   0.4199  
##  Residual                 0.6704   0.8188  
## Number of obs: 396, groups:  condition:id, 121; id, 65
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)          1.38757    0.09252  99.37271  14.998  < 2e-16 ***
## conditiondrug        0.67315    0.11522  56.46129   5.842 2.68e-07 ***
## time                 0.10055    0.04915 287.77919   2.045   0.0417 *  
## conditiondrug:time  -0.11214    0.07481 291.39762  -1.499   0.1350    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.539              
## time         0.078 -0.053       
## cndtndrg:tm -0.051  0.058 -0.656
## [1] "###Neurocaps MCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time +  (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 428.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8776 -0.6184 -0.1551  0.2547  3.1725 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.08263  0.2875  
##  id           (Intercept) 0.14869  0.3856  
##  Residual                 0.71649  0.8465  
## Number of obs: 156, groups:  condition:id, 48; id, 26
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)          1.15054    0.12925  30.42153   8.902 5.65e-10 ***
## conditiondrug        0.83342    0.16704  20.31452   4.989 6.74e-05 ***
## time                -0.05272    0.07641 102.21618  -0.690 0.491729    
## conditiondrug:time  -0.48453    0.12643 112.93184  -3.832 0.000209 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.498              
## time         0.029 -0.018       
## cndtndrg:tm -0.021 -0.043 -0.607
## [1] "###Neurocaps OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time +  (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 346
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9132 -0.4724 -0.1693  0.3446  3.3533 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.1120   0.3347  
##  id           (Intercept) 0.3614   0.6012  
##  Residual                 0.4676   0.6838  
## Number of obs: 139, groups:  condition:id, 40; id, 20
## 
## Fixed effects:
##                    Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)         1.52428    0.17328 25.16252   8.797 3.77e-09 ***
## conditiondrug       0.26377    0.16006 15.67478   1.648   0.1193    
## time                0.03842    0.07248 93.95191   0.530   0.5973    
## conditiondrug:time -0.25231    0.10443 97.08704  -2.416   0.0176 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.428              
## time         0.001  0.000       
## cndtndrg:tm  0.002  0.026 -0.693
## [1] "###tACS OCR###"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: paste("value ~ condition * time +  (1|id/condition)")
##    Data: long_data
## 
## REML criterion at convergence: 408.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6896 -0.5276 -0.1386  0.2052  4.5287 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  condition:id (Intercept) 0.02208  0.1486  
##  id           (Intercept) 0.42820  0.6544  
##  Residual                 0.49879  0.7063  
## Number of obs: 165, groups:  condition:id, 44; id, 22
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)          1.27609    0.16272  25.57989   7.842 2.87e-08 ***
## conditiondrug        0.39445    0.11926  15.04231   3.307  0.00477 ** 
## time                -0.04465    0.06858 111.57509  -0.651  0.51634    
## conditiondrug:time  -0.09104    0.09924 114.56134  -0.917  0.36086    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndtnd time  
## conditindrg -0.360              
## time         0.008 -0.007       
## cndtndrg:tm -0.005  0.016 -0.689
extract_slopes_beh <- function(measure, this_label, idps, prefix = 'dcr_'){
  #same as above, but for behavioral data
  #prefix will be dcr_tpre_ for tdcs/tacs and dcr_ for neurocaps 
  #measure will be box_rt, craving_rt, or craving
  #measure <- 'craving'
  #prefix <- 'dcr_tpre_'
  #this_label <- 'Craving'
  #idps <- idps_tdcs
  
  
  #idps_tacs$dcr_tpre_response_craving_[0-7]
  #idp_tdcs
  #dcr_tpre_response_box_rt_0
  #dcr_tpre_response_craving_rt_0
  #idps_neurocaps_ocr$dcr_response_box_craving[0-7]

  column_suffixes <- 0:7
  
  these_cols <- c('id', paste0(prefix, 'response_', measure, '_', column_suffixes))
  library(reshape2)
  one_dataset <- idps[, these_cols]
  long_data <- melt(one_dataset, id.vars = c('id'))
  long_data$variable <- as.character(long_data$variable)
  long_data$condition <- NA
  long_data$number <- substr(long_data$variable, nchar(long_data$variable), nchar(long_data$variable))
  long_data$condition[long_data$number %in% c('0', '2', '4', '6')] <- 'neutral'
  long_data$condition[long_data$number %in% c('1', '3', '5', '7')] <- 'drug'
  #put neutral in the intercept
  long_data$condition <- factor(long_data$condition, levels = c('neutral', 'drug'))
  #time = block number, just like for the imaging variables
  long_data$time <- NA
  long_data$time[long_data$number %in% c('0', '1')] <- 1
  long_data$time[long_data$number %in% c('2', '3')] <- 2
  long_data$time[long_data$number %in% c('4', '5')] <- 3
  long_data$time[long_data$number %in% c('6', '7')] <- 4

  #mean center on time
  long_data$time <- long_data$time - mean(long_data$time)

  library(lme4)
  library(lmerTest)
  library(ggplot2)
  library(sjstats) #for icc of mixed effects models
  
  
  
  all_results <- NULL
  for (this_id in unique(long_data$id)){  
    #this_id <- 'AR316'
    this_sub_data <- long_data[long_data$id == this_id,]
    neutral_string <- paste0(this_label, '_slopeNeutral')
    meth_string <- paste0(this_label, '_slopeMeth')
    delta_string <- paste0(this_label, '_slopeMethMinusNeutral')
    if (sum(!(is.na(this_sub_data$value[this_sub_data$condition == 'neutral']))) > 1 &
        sum(!(is.na(this_sub_data$value[this_sub_data$condition == 'drug']))) > 1){
      #subject made at least two ratings per condition
      this_lm <- lm(paste('value ~ condition * time'), data = this_sub_data)
      
      these_results <- data.frame(list(id = this_id, neutral_string = this_lm$coefficients['time'],
                                meth_string = this_lm$coefficients['time'] + this_lm$coefficients['conditiondrug:time'],
                                delta_string = this_lm$coefficients['conditiondrug:time']))
      names(these_results) <- c('id', neutral_string, meth_string, delta_string)
      
    } else{
      #not enough ratings to estimate two slopes
      these_results <- data.frame(list(id = this_id, neutral_string = NA,
                                meth_string = NA,
                                delta_string = NA))
      names(these_results) <- c('id', neutral_string, meth_string, delta_string)

        }
    all_results <- rbind(all_results, these_results)
  }
  

  return(all_results)
   
}

Corrplot including individual slopes and self-report

extract_one_slopes <- function( roi = '1', this_label = 'VMPFC', prefix = 'dcr_tpre_', idps = idps_tdcs){
  column_prefixes <- c('stats_tdcsprelim_drug.r11.0.coef_mean_',
           'stats_tdcsprelim_drug.r12.0.coef_mean_',
           'stats_tdcsprelim_drug.r13.0.coef_mean_',
           'stats_tdcsprelim_drug.r14.0.coef_mean_',
           'stats_tdcsprelim_neutral.r11.0.coef_mean_',
           'stats_tdcsprelim_neutral.r12.0.coef_mean_',
           'stats_tdcsprelim_neutral.r13.0.coef_mean_',
           'stats_tdcsprelim_neutral.r14.0.coef_mean_')
  this_roi <- c('id', 'motion', paste0(prefix, column_prefixes, roi))
  one_dataset <- idps[, this_roi]
  long_data <- melt(one_dataset, id.vars = c('id', 'motion'))
  long_data$condition <- NA
  long_data$condition[grepl('neutral', long_data$variable)] <- 'neutral'
  long_data$condition[grepl('drug', long_data$variable)] <- 'drug'
  #put neutral in the intercept
  long_data$condition <- factor(long_data$condition, levels = c('neutral', 'drug'))
  long_data$time <- NA
  long_data$time[grepl('r11', long_data$variable)] <- 1
  long_data$time[grepl('r12', long_data$variable)] <- 2
  long_data$time[grepl('r13', long_data$variable)] <- 3
  long_data$time[grepl('r14', long_data$variable)] <- 4
  #mean center on time
  long_data$time <- long_data$time - mean(long_data$time)
  
  #this_lme <- lmer(paste('value ~ condition * time + (1|id)'), data = long_data)
  #for checking for NA's, looks like it's all good now
  #print(long_data)
  all_results <- NULL
  for (this_id in unique(long_data$id)){  
    #this_id <- 'AR316'
    this_lm <- lm(paste('value ~ condition * time '), data = long_data[long_data$id == this_id,])
    neutral_string <- paste0(this_label, '_slopeNeutral')
    meth_string <- paste0(this_label, '_slopeMeth')
    delta_string <- paste0(this_label, '_slopeMethMinusNeutral')
    these_results <- data.frame(list(id = this_id, neutral_string = this_lm$coefficients['time'],
                                meth_string = this_lm$coefficients['time'] + this_lm$coefficients['conditiondrug:time'],
                                delta_string = this_lm$coefficients['conditiondrug:time']))
    names(these_results) <- c('id', neutral_string, meth_string, delta_string)
    all_results <- rbind(all_results, these_results)
  }
  
  

  return(all_results)
}


all_slopes <- extract_one_slopes('1', 'VMPFC', prefix = 'dcr_tpre_', idps = idps_tdcs)
all_slopes <- merge(all_slopes, extract_one_slopes('3', 'LSTG', prefix = 'dcr_tpre_', idps = idps_tdcs))
all_slopes <- merge(all_slopes, extract_one_slopes('4', 'RSTG', prefix = 'dcr_tpre_', idps = idps_tdcs))
all_slopes <- merge(all_slopes, extract_one_slopes('5', 'LVSTRI', prefix = 'dcr_tpre_', idps = idps_tdcs))
all_slopes <- merge(all_slopes, extract_one_slopes('6', 'RVSTRI', prefix = 'dcr_tpre_', idps = idps_tdcs))
all_slopes <- merge(all_slopes, extract_one_slopes('8', 'RAmy', prefix = 'dcr_tpre_', idps = idps_tdcs))






vas_data <- read.csv('../paper-dcr-temporaldynamics/data/MethVASData.csv')
vas_data <- vas_data[, c('record_id', 'redcap_event_name', 'mcs_vas', 'mcs_vas_2')]

names(vas_data) <- c('id', 'visit', 'craving', 'control')

library(reshape2)
vas_data_1 <- vas_data[vas_data$visit == 'before_pre_fmri_arm_1', c('id', 'craving', 'control')]
vas_data_2 <- vas_data[vas_data$visit == 'after_pre_fmri_arm_1', c('id', 'craving', 'control')]

vas_data_wide <- merge(vas_data_1, vas_data_2, by = 'id')

names(vas_data_wide) <- c('id', 'craving_pre', 'control_pre', 'craving_post', 'control_post')
vas_data_wide$craving_delta <- vas_data_wide$craving_post - vas_data_wide$craving_pre
vas_data_wide$control_delta <- vas_data_wide$control_post - vas_data_wide$control_pre

merged_data <- merge(all_slopes, vas_data_wide, all.x = TRUE)



last_use_data <- read.csv('../paper-dcr-temporaldynamics/data/ADUQ_2019-06-12_1407.csv')
merged_data <- merge(merged_data, last_use_data[, c('id', 'aduq_20a')], all.x = TRUE)
names(merged_data)[names(merged_data) == 'aduq_20a'] <- 'DaysSinceLastUse'




tableone_data <- read.csv('../paper-dcr-temporaldynamics/data/Table1Database66_ver2.csv')
merged_data <- merge(merged_data, tableone_data, all.x = TRUE)

idps_tdcs$methvsneutral_craving_selfreport_insidescanner <- (idps_tdcs$dcr_tpre_response_craving_1 + idps_tdcs$dcr_tpre_response_craving_3 +
  idps_tdcs$dcr_tpre_response_craving_5 + idps_tdcs$dcr_tpre_response_craving_7 - 
    idps_tdcs$dcr_tpre_response_craving_0 - idps_tdcs$dcr_tpre_response_craving_2 - idps_tdcs$dcr_tpre_response_craving_4 - idps_tdcs$dcr_tpre_response_craving_6) / 4


merged_data <- merge(merged_data, extract_slopes_beh('craving', 'Craving', idps_tdcs, prefix = 'dcr_tpre_'), all.x = TRUE)

merged_data <- merge(merged_data, idps_tdcs[, c('id', 'methvsneutral_craving_selfreport_insidescanner')])
names(merged_data)
##  [1] "id"                                                "VMPFC_slopeNeutral"                                "VMPFC_slopeMeth"                                   "VMPFC_slopeMethMinusNeutral"                       "LSTG_slopeNeutral"                                
##  [6] "LSTG_slopeMeth"                                    "LSTG_slopeMethMinusNeutral"                        "RSTG_slopeNeutral"                                 "RSTG_slopeMeth"                                    "RSTG_slopeMethMinusNeutral"                       
## [11] "LVSTRI_slopeNeutral"                               "LVSTRI_slopeMeth"                                  "LVSTRI_slopeMethMinusNeutral"                      "RVSTRI_slopeNeutral"                               "RVSTRI_slopeMeth"                                 
## [16] "RVSTRI_slopeMethMinusNeutral"                      "RAmy_slopeNeutral"                                 "RAmy_slopeMeth"                                    "RAmy_slopeMethMinusNeutral"                        "craving_pre"                                      
## [21] "control_pre"                                       "craving_post"                                      "control_post"                                      "craving_delta"                                     "control_delta"                                    
## [26] "DaysSinceLastUse"                                  "Age"                                               "Education"                                         "BMI"                                               "Age.of.Meth.use.onset..years.old."                
## [31] "Duration.of.Meth.use.at.least.once.a.week..years." "Cost.of.Meth..dollar.per.month."                   "Dose.of.Meth..gram.per.day."                       "History.of.Meth.injection.b..n....."               "Days.of.Meth"                                     
## [36] "Days.of.Alcohol"                                   "Days.of.Alcohol.intoxication"                      "Days.of.Heroin"                                    "Days.of.Methadone"                                 "Days.of.Other.opioids"                            
## [41] "Days.of.Barbiturate"                               "Days.of.Sedative"                                  "Days.of.Cocaine"                                   "Days.of.Cannabis"                                  "Days.of.Hallucinogens"                            
## [46] "Days.of.Inhalants"                                 "Duration.of.current.abstinence..days."             "Meth.Cue.Reactivity.Screening.score..0.100."       "Meth.Withdrawal.Scale.score..0.24."                "Motion..Rest.fMRI..Pre.tDCS"                      
## [51] "Motion..Task.fMRI..Pre.tDCS"                       "Motion..Rest.fMRI..Post.tDCS"                      "Motion..Task.fMRI..Post.tDCS"                      "Craving_slopeNeutral"                              "Craving_slopeMeth"                                
## [56] "Craving_slopeMethMinusNeutral"                     "methvsneutral_craving_selfreport_insidescanner"
names(merged_data)[names(merged_data) == 'Meth.Cue.Reactivity.Screening.score..0.100.'] <- 'BaselineCueReactivity'
names(merged_data)[names(merged_data) == 'Duration.of.Meth.use.at.least.once.a.week..years.'] <- 'MethUseDuration'
names(merged_data)[names(merged_data) == 'Cost.of.Meth..dollar.per.month.'] <- 'MethCost'

names(merged_data)[names(merged_data) == 'craving_post'] <- 'Craving_post'
names(merged_data)[names(merged_data) == 'craving_pre'] <- 'Craving_pre'
names(merged_data)[names(merged_data) == 'craving_delta'] <- 'Craving_delta'

to_plot <- c('Age', 'MethUseDuration', 'MethCost',"DaysSinceLastUse",
            'BaselineCueReactivity',  "Craving_pre", "Craving_post", "Craving_delta",
            "VMPFC_slopeNeutral", "VMPFC_slopeMeth",  "VMPFC_slopeMethMinusNeutral",
            "LSTG_slopeNeutral", "LSTG_slopeMeth","LSTG_slopeMethMinusNeutral",
            "RSTG_slopeNeutral", "RSTG_slopeMeth", "RSTG_slopeMethMinusNeutral",
            "LVSTRI_slopeNeutral", "LVSTRI_slopeMeth", "LVSTRI_slopeMethMinusNeutral",
            "RVSTRI_slopeNeutral", "RVSTRI_slopeMeth", "RVSTRI_slopeMethMinusNeutral",
            "RAmy_slopeNeutral", "RAmy_slopeMeth", "RAmy_slopeMethMinusNeutral")

to_plot <- c('Age', 'MethUseDuration', 'MethCost',"DaysSinceLastUse",
            "Craving_pre", "Craving_post", "Craving_delta", 
            'methvsneutral_craving_selfreport_insidescanner',
            'Craving_slopeMethMinusNeutral',
            "VMPFC_slopeMethMinusNeutral",
            "LSTG_slopeMethMinusNeutral",
            "RSTG_slopeMethMinusNeutral",
            "LVSTRI_slopeMethMinusNeutral",
            "RVSTRI_slopeMethMinusNeutral",
            "RAmy_slopeMethMinusNeutral")

library(corrplot)
## corrplot 0.84 loaded
this_matrix <- cor(merged_data[, to_plot], use = 'pairwise.complete.obs')
pvals_raw <- cor.mtest(merged_data[, to_plot])

#get FDR corrected p-values to use for plotting
pvals_fdr <- pvals_raw$p
pvals_fdr[upper.tri(pvals_fdr)] <- p.adjust(pvals_raw$p[upper.tri(pvals_raw$p)])
pvals_fdr[lower.tri(pvals_fdr)] <- p.adjust(pvals_raw$p[lower.tri(pvals_raw$p)])


col2 <- colorRampPalette(rev(c("#67001F", "#B2182B", "#D6604D", "#F4A582",
                           "#FDDBC7", "#FFFFFF", "#D1E5F0", "#92C5DE",
                           "#4393C3", "#2166AC", "#053061")))


png('corrplot_slopes_sxs.png', width = 1000, height = 1000)
corrplot.mixed(this_matrix, upper = 'ellipse', lower = 'number', tl.pos = 'lt', tl.cex = 1,#, lower.col = "black",
  upper.col = col2(50), lower.col = col2(50), p.mat = pvals_fdr, sig.level = 0.05)
dev.off()
## png 
##   2